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SUMMARY 


This  report  gives  the  results  of  a  Monte  Carlo  calculation  of  the 
energy  spectrum  and  angular  distribution  of  gamma  radiation  emitted  by  a 
collimated  cobalt-60  source.  It  is  an  extension  of  a  previous  report 
from  the  University  of  Illinois  (K.  Preiss,  "Monte  Carlo  Calculation  of 
Self  Shielding  by  Encapsulated  Gamma  Ray  Sources,"  Univ.  of  Ill.  Rpt. 
NRSS-3,  August  1966). 

The  program  developed  is  flexible  enough  to  give  results  for  a  wide 
variety  of  geometric  situations  of  source  container  and  collimator;  how¬ 
ever,  the  results  of  the  specific  calculation  are  for  the  particular  one 
used  by  the  University  of  Illinois  in  its  shielding  experimentation. 

Minor  simplifications  ir,  the  precise  geometry  have  been  made  for  ease  of 
programming,  but  such  simplifications  are  not  of  a  nature  to  change  the 
results  to  any  appreciable  extent. 

In  the  Monte  Carlo  calculation,  the  major  simplification  in  technique 
involves  the  use  of  a  fictitious  interaction  to  facilitate  the  handling 
of  some  of  the  boundary  problems.  In  essence,  if  this  interaction  is 
properly  introduced,  two  media  having  tabulated  attenuation  coefficients 
which  are  slightly  different  can  be  made  to  have  the  same  total  interac¬ 
tion  coefficient  for  the  purpose  of  calculating  path  lengths,  anc  the 
points  of  boundary  crossing  between  these  two  media  need  not  be  calculated 
(This  technique  seems  to  have  been  first  mentioned  by  Woodcock  et  al. ,  U.S 
AEC  Rpt.  TID  4500,  May  1965. ) 

The  output  data  are  given  in  terms  of  13  energy  groups,  from  1.4  MeV 
down  to  0.1  MeV,  and  3  angle  groups.  The  angle  groups  involve:  first,  a 
narrow  cone  in  which  all  radiation  passing  through  the  collimator  without 
interaction  must  be  included;  a  somewhat  wider  cone  which  is  established 
by  a  secondary  coll imator  placed  on  top  the  first  tc  minimize  lip  scatter 
problems;  and  the  remainder  which  includes  all  lip  scatter  out  to  90°  from 
the  axis  of  the  collimator.  'he  first  angular  group  contains  about  80% 
of  the  emergent  photons  and  the  other  two  groups  contain  about  10%  each. 

For  the  spectrum  in  the  first  angular  interval,  scattered  photons 
constitute  16%.  The  low  energy  spectrum  rises  to  a  peak  in  the  energy 


/•*" 


interval  0.2  to  0.3  MeV.  The  peak  results  primarily  from  the  fact  that  the 
single-scattering  cutoff  energies  for  bolh  source  energies  lie  in  this  in¬ 
terval.  The  spectrum  results  almost  entirely  from  scatterings  in  the 
source  and  its  capsule,  with  very  little  contribution  from  collimator 
scatterings. 

The  results  are  in  qualitative  agreement  with  previous  calculations 
published  by  Cormack  and  Johns  for  a  somewhat  similar  situation  (British 
Journal  of  Radiology,  31,  No.  369,  497,  1958).  Agreement  is  also  good 
with  experimental  measurements  of  Costrell  on  encapsulated  sources 
(Health  Physics,  8,  261,  1962),  except  at  very  low  energies. 

Modifications  of  the  emergent  spectrum  by  introduction  of  lead  filter 
plates  across  the  collimator  are  included. 
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ABSTRACT 


A  Monte  Carlo  calculation  of  the  energy  spectrum  emitted  by  a  col¬ 
limated,  Co-60  source  is  reported.  The  collimator  assumed  is  very 
similar  to  one  currently  being  used  in  experimental  shielding  studies 
at  the  University  of  Illinois.  Radiation  emerging  from  the  collimator 
is  classified  into  13  equal-length  energy  intervals  between  0.1  and 
1.4  MeV  ,  and  into  3  direction  intervals  for  angles  relative  to  the 
collimator  axis  between  0°  and  90°  .  Comparisons  are  made  w; th 

experimental  and  theoretical  work  reported  by  others. 
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I,  INTRODUCTION 

This  report  gives  the  results  of  a  Monte  Carlo  calculation  of  the 

energy  spectrum  and  angular  distribution  of  gamma  radiation  emitted  by 

a  collimated  Co-60  source.  It  is  an  extension  of  earlier  work  reported 

(1) 

by  Preiss  for  uncollimated,  encapsulated  sources. 

Preiss  cites  several  instances  where  discrepancies  between  theo¬ 
retical  calculations  and  corresponding  experimental  results  have  been 
traced  to  neglect  of  source  spectrum  degradation  due  to  scattering 
either  in  the  source  material  or  in  material  near  the  source.  The 
spectra  in  this  report  are  to  be  used  as  input  for  calculations  of 
gamma-ray  penetration  through  ribbed  slabs.  The  results  of  these  cal¬ 
culations  in  turn  will  be  compared  with  results  of  experimental  studies 
of  ribbed-slab  shielding  properties  being  conducted  at  the  University  of 
Illinois. 

The  collimated  source  for  which  the  present  calculation  was  made  is 
shown  in  Figure  1.  The  source  being  used  in  the  experimental  work  men¬ 
tioned  above  differs  from  the  collimated  source  in  Figure  1  in  that  the 
experimental  source  capsule  is  held  in  place  by  a  collar  made  of  Firth 
Heavy  Metal.  This  alloy  consists  of  about  90%  tungsten  with  the  remain 
ing  10%  consisting  of  copper  and  nickel.  In  Figure  1,  the  collar  is  re 
garded  as  "black."  The  approximation  thus  made  is  not  expected  to  have  a 
significant  influence  on  the  usefulness  of  the  calculated  results. 

The  Monte  Carlo  calculation  on  which  the  data  in  this  report  are 
based  is  described  in  Section  II.  The  mass  attenuation  coefficients  used 
in  the  calculation  are  given  in  Appendix  B  along  with  a  description  of 
the  procedure  employed  to  evaluate  interaction  coefficients  for  cobalt. 
The  computer  program  for  the  Monte  Carlo  calculation  is  described  in 
Appendix  C  . 

Some  experimental  measurements  of  spectra  from  collimated  Co-60 

sources  have  been  reported. Among  these,  one  of  the  sources 

(2) 

considered  by  Costrell  corresponds  fairly  closely  to  the  collimated 
source  in  Figure  1  ,  We  compare  with  one  of  his  spectra  in  Section  III. 
Also  in  Section  III,  our  results  are  compared  with  the  Monte  Carlo  work 
of  Preiss^  and  with  a  calculation  by  Cormack  and  Johns.  ^ 
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Frequently,  in  order  to  reduce  the  importance  of  low  energy  scattered 
radiation  relative  to  primary  radiation,  the  source  spectrum  is  filtered 
by  placing  a  thin  sheet  of  lead  between  the  source  and  the  target.  In 
Appendix  A,  spectra  which  have  been  filtered  in  this  way  are  presented 
for  lead  thicknesses  of  1/8  and  1/4  inch. 

II.  DESCRIPTION  OF  THE  CALCULAT I ONAL  BASIS  OF  THE  MONTE  CARLO  PROGRAM 

The  axis  of  the  collimator  was  taken  as  the  Z-a-<is  of  a  cartesian 
coordinate  system.  Photon  positions  were  then  given  in  terms  of  rectangu¬ 
lar  coordinates  and  photon  directions  were  specified  by  direction  cosines 
relative  to  each  of  the  coordinate  axes. 

The  Co-60  was  assumed  to  be  uniformly  distributed  in  the  cobalt 
metal  (Figure  1).  At  the  beginning  of  the  calculation,  a  finite  set  of 
source  points  which  were  assumed  to  represent  the  uniform  distribution 
were  selected.  The  selection  was  made  by  dividing  the  source  region  into 
a  number  of  equal  volume  increments  and  picking  a  single  source  point  in 
each  increment.  Because  of  cylindrical  symmetry,  the  volume  increments 
were  chosen  in  the  shape  of  annuli.  They  were  defined  by  dividing  the 
source  region  into  a  number  of  equal  thickness  layers  all  parallel  to 
the  plane  z  =  0  ,  and  then  further  dividing  each  layer  into  several 
concentric  rings.  The  source  point  for  a  particular  volume  increment 
was  selected  so  that  it  was  located  on  the  XZ-plane  with  x  >  0  ,  half¬ 
way  between  the  top  and  bottom  planes,  and  halfway  between  the  inner 
and  outei  radii  of  the  volume.  As  the  calculation  progressed,  source 
points  were  then  systematically  selecttd  from  this  predetermined  list 
so  that  for  each  of  the  source  energies  1.17  and  1.33  MeV  ,  the  same 
number  of  photon  histories  originated  at  each  source  point. 

Initial  photon  directions  were  sampled  randomly  from  an  isotropic 
distribution. 

The  basic  interaction  considered  in  the  calculation  other  than 
absorption  was  Compton  scattering.  Absorption  was  taken  primarily  as 
photoelectric;  but  for  the  short  energy  interval  1.022  to  1.33  MeV  , 
the  pair  production  interaction  was  also  treated  as  absorptive. 

In  addition,  we  followed  the  procedure  first  reported  by  Woodcock, 
et  al.  ^  and  subsequently  used  by  Preiss, ^  of  introducing  a  fictitious 
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interaction  to  facilitate  the  handling  of  some  of  the  boundary  problems. 
This  procedure  may  be  described  as  follows:  Consider  two  media  with  an 
interface  between  them.  Let  p^  and  p^  be  the  linear  attenuation 
coefficients  for  medium  1  and  medium  2  respectively,  further,  suppose 
p  >  p  .  If  a  photon  interacts  in  medium  2  and  is  not  absorbed,  then 
we  sample  a  path  length  to  the  next  interaction  point  using  the  equation, 


where  s  is  the  path  length  and  p  is  a  random  number  between  0  and 
1  .  If  the  path  carries  the  photon  across  the  interface  into  medium  1, 
we  write, 

^2  '  ^2  ”  ^  T  ^al  4  ^sl  ’ 

where 

**1  =  »*al  +  wsl  ’ 

and  where  p  and  p  are  the  linear  absorption  and  scattering  coef- 
al  si 

ficients  for  medium  1.  Three  possible  interactions  are  then  considered 
instead  of  the  usual  two.  The  first,  with  probability 

1^2  ~ 

is  a  fictitious  interaction  in  which  there  is  no  energy  loss,  no  change 
of  direction,  and  no  aosorption.  The  second,  with  probability 


is  a  Compton  scattering  interaction.  The  third,  with  probability 


p<2 


is  a  photoelectric  absorption.  Ihe  introduction  of  the  fictitious 
interaction  in  medium  1  allows  us  to  use  the  linear  attenuation  coef¬ 
ficient  p^  for  sampling  path  lengths  in  both  medium  1  and  medium  2. 
Thus,  it  is  no  longer  necessary  to  compute  the  coordinates  of  the 
crossing  i-oint  every  time  the  photon  crosses  an  interface  or  to  check 
which  of  the  various  bounding  surfaces  has  been  crossed  first. 

In  the  present  calculation,  the  above  orocedure  was  used  in  the 
source  and  capsule,  and  separately  in  the  collimator.  In  the  source 
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and  capsule,  cobalt  was  medium  2  and  iron  was  medium  1.  In  the  collimator, 
lead  was  medium  2  and  iron  was  medium  1.  (In  principle,  it  would  have 
been  possible  to  base  transport  in  the  whole  assembly  on  lead;  however, 
this  would  have  resulted  in  a  very  large  number  of  fictitious  interactions 
in  the.  source  and  capsule,  with  attendant  increase  in  computer  time.) 

The  air  column  along  the  axis  of  the  collimator  was  treated  as  a 
vacuum.  Thus,  when  a  photon  entered  this  region,  it  was  moved  along 
its  trajectory  until  it  either  re-entered  the  collimator  wall  or  source 
capsule,  or  until  it  passed  out  the  end  of  the  collimator.  For  this 
reason,  it  was  necessary  to  compute  the  coordinates  of  the  crossing 
point  when  the  photon  crossed  one  of  the  surfaces  bounding  the  air 
column. 


The  calculation  of  crossing  point  coordinates  is  very  simple  when 


the  photon  crosses  a  plane  z  =  z 


If  the  photon  starts  at  the  point 


(x, ,  y. ,  z.)  and  moves  along  a  path  with  direction  cosines  u  ,  u  ,  ar 

ill  or-  x  y 

uz,  then  the  distance  s  which  the  photon  travels  to  reach  the  plane 
z  -  %  18  ■ 


The  coordinates  of  the  crossing  point  are  given  by, 


X  = 

xi 

+ 

su  , 
X 

y  = 

yi 

+ 

su  , 

y 

z  = 

Zb 

• 

The  calculation  of  crossing  point  coordinates  is  more  complicated 
if  the  photon  leaves  the  point  (x^,  y  ,  z^)  and  crosses  a  cylindrical 
boundary  with  radius  r  .  If  x  and  y  are  the  x-  and  y-  coordi¬ 
nates  of  the  crossing  point,  then 

2  2  2 
r  =  x  +  y 

Hut  x  and  y  are  given  in  terms  of  x  and  y  by  the  equation, 

x  =  x,  +  su  , 
i  x 

y  =  yt  +  suy  - 

where  s  is  the  distance  from  the  point  (x^,  y  ,  z^)  to  the  crossing 
point.  Therefore,  s  must  satisfy  the  equation, 
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2  2  2 
r  =  (x  r  su  )“  (y  -  su  ) 
ix  i  y 

2  2 

=  x  “  *  y  r  2s ( u  x  *uy) 
i  i  xi  y  i 

2  2  2 

+  s  ( u  -  u  ) 
x  y 


Solving  this  equation,  we  iind. 


/ 2  2  2,  2  2 

~(ux  +  u  y  )  ±  \/(ux  -  u  v  )  Tfu  -u  )  (r  ~  x  -  y.  ) 

x  i  v  l  V  x  i  V  l  x  v  l  l 


,  2  2, 
(u  *  U  ) 
x  y 


The  sign  of  the  square  root  must  be  chosen  according  to  the  rules: 


U  X.  + 
X  1 

U  V. 

y  i 

2 

0. 

2 

r 

> 

2 

x  + 

l 

2 

yi 

plus 

sign; 

U  X. 

X  1 

u  y 

y  i 

< 

0, 

o 

r 

2 

2 

X 

i 

2 

yi 

plus 

sign; 

U  X.  * 
X  1 

u  y. 
y  v 

< 

0 

2 

r 

< 

2 

x  + 

l 

2 

yi 

minus 

sign. 

Once  s  ha»  been  calculated,  the  coordinates  of  the  crossing  point  are, 

X  =  x.  u  s 
1  x 

y  =  y.  -r  u  s  , 
l  y 

Z  =  Z  t  u  s  . 

1  Z 

Photon  histories  were  terminated  under  the  following  conditions: 

(1)  if  a  photon  entered  a  black  region;  (2)  if  at  a  given  interaction 
point  a  photoelectric  interaction  was  selected,  (2)  if  the  photon  energy 
dropped  below  0.1  MeV  in  iron  or  cobalt,  or  below  0.2  MeV  in  lead; 

or  (4)  if  the  photon  left  the  collimator. 

When  a  photon  emerged  from  the  end  of  the  collimator,  eithei  from 

the  air  column  or  from  the  surrounding  iron  and  lead  a  score  was  recorded. 

The  direction  cosine  u of  the  emergent  photon  was  classified  into  one 

of  the  intervals  0.996  s.  u  ^  1.0  0.96i  u  <  0,996.  or  0  S  u  <  0,961 

z  z  z 

and  its  energy  E  was  classified  into  one  of  13.  0. l"MeV  energy  inter¬ 
vals  between  0.1  and  1.4  MeV  .  The  scoring  procedure  consisted  of 
adding  the  statistical  weight  of  the  photon  (always  unity  in  the  present 
calculation)  to  a  r  nning  total  for  the  appropriate  energy  angular 
interval. 


The  three  angular  groups  defined  above  were  of  interest  with  respect 

to  the  University  of  Illinois  24-curie  (nominal)  experimental  source. 

All  uncollided  photons  which  emerged  without  penetrating  any  part  of  the 

lip  of  the  collimator  had  to  have  direction  cosines  u^  >  0.996,  while 

photons  which  emerged  with  u  <  0.996  had  to  penetrate  some  portion  of 

z 

the  lip.  For  the  experimental  source,  most  of  the  photons  which  emerged 
with  u^  <  0.961  had  to  penetrate  an  auxiliary  lead  collimator  which  was 
placed  above  the  main  collimator,  as  indicated  in  Figure  1. 

III.  RESULTS  AND  COMPARISON  WITH  OTHER  WORK 

The  results  of  the  calculation  described  in  the  preceding  section 
are  given  in  Table  1.  The  data  are  based  on  one  million  histories. 

This  large  number  of  histories  was  required  because  of  the  small  angle 
of  collimation  and  because  of  the  fairly  large  number  (39)  of  energy- 
angle  groups  used  in  the  calculation.  About  63  minutes  on  the  IBM  7094 
at  the  University  of  Illinois  were  required  for  the  calculation. 

An  analytical  calculation  of  the  number  of  uncollided  photons  which 
one  would  expect  to  emerge  from  the  collimator  was  made.  Penetration 
of  the  lip  of  the  collimator  was  neglected.  Using  the  normalization  of 
Table  1,  the  analytical  calculation  predicted  781  photons.  From 
Table  1,  we  see  that  948  photons  emerge  in  the  energy  intervals  1.1 
to  1.2  MeV  and  1.3  to  1.4  MeV  and  the  angular  interval  0.966  £  u 

z 

£1.0  .  These  photons  may  be  regarded  as  uncollided.  (An  approximate 

estimate  of  the  number  of  scattered  photons  contributing  to  the  above- 
mentioned  energy  intervals  appears  to  be  between  10  and  20  ,)  This 
means  that  approximately  17%  of  all  uncollided  photons  emerging  from 
the  collimator  have  penetrated  some  portion  of  the  lip. 

The  spectrum  of  gamma  radiation  emerging  from  the  collimator  in 
the  angular  interval  0.996  £  u  £  1.0  is  shown  in  Figure  2.  Scattered 
photons  contribute  16%  of  the  area  under  the  histogram.  The  low  energy 
spectrum  rises  to  a  peak  in  the  energy  interval  0.2  to  0.3  MeV  . 

The  peak  results  primarily  from  the  fact  that  the  single-scattering  cutoff 
energies  for  both  source  energies  lie  in  this  interval.  The  spectrum 
decreases  rapidly  below  0.2  MeV  because  of  the  increasing  importance 
of  the  photoelectric  effect  at  low  photon  energies. 
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TABLE  1 


Number  of  photons  p/>r  0. 1  MeV  per  500,000  disintegrations 
of  Co-60.  The  standard  deviation  for  an  entry  in  the 
table  can  be  estimated  by  computing  the  square  root. 


\  Angular 

Interval 
Energy  \ 
Interval 
(MeV)  \ 

0.996  £  u  £1.0 

2 

0.961  £  u  <  0.996 

2 

0  £  u  <0. 961 
z 

493 

7 

0 

3 

18 

l 

1.2  >  E  £  1.1 

455 

31 

8 

1.1  >  E  £  1.0 

13 

36 

26 

1.0  >  E  £  0.9 

6 

1 

17 

19 

0.9  >  E  a  0.8 

ii 

8 

25 

0.8  >  E  £  0.7 

14 

4 

19 

0. 7  >  £  £  0.6 

12 

9 

14 

0.6  >  E  £  0.5 

19 

8 

14 

0.5  >  E  £  0.4 

22 

3 

14 

0.4  >  E  £0.3 

30 

1 

11 

0.3  >  E  £  0.2 

36 

5 

4 

0.2  >  E  £  0.1 

16 

5 

3 

Sum 

1130 

152 

158 
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The  spectrum  in  Figure  2  results  almost  entirely  from  scatterings 
in  the  source  and  capsule.  This  can  be  seen  if  we  compare  with  one  of 
the  spectra  computed  by  Preiss^  for  an  uncollimated,  encapsulated 
source.  In  Figure  3,  Preiss's  spectrum  for  his  source  B  is  compared 
with  the  spectrum  in  Figure  2.  Preiss's  source  consisted  of  a  cobalt 
metal  disk  inside  a  steel  capsule.  The  cobalt  disk  hud  the  same  dimen¬ 
sions  as  the  disk  in  Figure  1.  The  capsule  had  the  same  radial  dimensions 
as  the  capsule  in  Figure  1,  but  its  length  was  only  0.5  cm,  Preiss's 

data  were  computed  for  0. 9  £  ur  S  1.0  ,  but  it  seems  reasonable  to 

z 

assume  that  the  shape  of  the  spectrum  was  practically  constant  over 

this  range  of  directions,  especially  in  view  of  the  similarity  of  Preiss's 

results  ^  for  0.9  ^  u  SI  and  0  S  u  <.0.1  ,  The  «ood  ugrseme"!” 

z  z 

between  our  spectrum  and  Preiss's  indicates  that  the  spectrum  in  Figure 
is  affectea  predominantly  by  the  top  half-centimeter  of  cc  e  e.tCTjsi  1  at e\I 
source. 

(2) 

Costrell  has  measured  spectra  for  a  variety  of  encapsulated, 

Co-60  sources,  both  collimated  and  uncollimated.  In  Figure  4,  the 
measured  spectrum  for  his  collimated  source  E2,5H  is  compared  with  the 
spectrum  presented  in  Figure  2.  His  collimator  was  made  of  lead,  was 
about  IS  inches  long,  and  had  a  diameter  varying  from  1  3/16  inches  at 
the  source  to  1  1/2  inches  at  the  mouth.  The  source  consisted  of  a 
disk  of  cobalt  metal  with  thickness  0.25  cm  and  diameter  1.7  cm.  It 
was  clad  in  aluminum  of  thickness  0.05  cm.  The  source  was  placed  in 
a  capsule  made  of  steel  and  tungsten.  The  thickness  of  the  source 
capsule  between  the  cobalt  and  the  end  of  the  collimator  was  0.035 
inches  of  steel.  The  E2.5H  source  was  selected  for  comparison  because 
this  thickness  was  almost  the  same  as  the  corresponding  thickness  for 
the  capsule  in  Figure  1,  Agreement  between  the  experimental  and  calcu¬ 
lated  results  is  quite  good  except  in  the  energy  interval  0.1  to 
0.2  MeV  . 

Figures  5  and  6  show  spectra  for  radiation  which  emerges  in  the 

angular  intervals  0.961  £  u  <  0,996  and  0  5  u  <  0,961,  Almost 

z  z 

all  the  radiation  emitted  in  these  angular  intervals  must  scatter  at 
least  once  from  the  walls  of  the  collimator.  (Uncollided  photons  may 
contribute  as  much  as  13%  of  the  histogram  area  in  Figure  5. )  For 
every  100  photons  which  emerge  from  the  collimator  uncollided,  i.e. 
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in  the  energy  intervals  1.1  to  1.2  MeV  and  1.3  to  1.4  MeV  and  in 
the  angular  interval  0.996  <.  ^  1.0  ,  about  16  photons  are  emitted 

in  each  of  the  direction  intervals  represented  in  Figures  5  and  6.  It 

is  interesting  to  note  that  most  of  these  photons  have  energies  greater 

than  0.5  MeV  . 

Cormack  and  Johns ^ ^  have  reported  a  single-scattering  calculation 
of  the  spectrum  of  radiation  scattered  from  the  walls  of  a  collimator. 
Their  collimator  was  made  of  lead  and  had  a  tapered  opening  with  a  mini¬ 
mum  radius  of  3  cm  near  the  source  and  a  maximum  radius  of  8  1/4  cm 
at  a  distance  of  60  cm  from  the  source.  They  calculated  the  spectrum 
at  a  point  on  the  axis  of  the  collimator  80  cm  from  the  source. 

In  Figure  7,  we  have  summed  the  spectra  presented  m  Figures  5  and  6 

and  normalized  to  one  uncollided  photon  emerging  from  the  collimator. 

The  shaded  portion  of  the  histogram  represents  the  contribution  from  the 
direction  interval  0.961  £  <  0.996  .  The  curve  is  the  spectrum 

calculated  by  Cormack  and  Johns.  It  rises  to  a  maximum  value  of  about 
2.5  at  an  energy  of  1.13  MeV  and  then  cuts  off  sharply.  In  the  Monte 
Carlo  calculation,  the  sharp  peak  is  smeared  out  over  several  energy 
intervals,  largely  because  of  use  of  a  finite  sized  source  and  detector. 

Because  of  considerable  geometric  differences,  we  expect  at  best 
only  qualitative  agreement  between  the  results  of  Cormack  and  Johns  and 
the  results  in  this  report.  Consequently,  the  near  equality  between  the 
area  under  the  curve  and  the  area  of  the  shaded  histogram  is  probably 
somewhat  fortuitous.  The  largest  angle  which  ladiation  directions  could 
have  with  respect  to  the  axis  of  the  collimator  and  be  included  in  the 
calculation  of  Cormack  and  Johns  was  about  22°  .  The  largest  angle 
which  was  included  in  the  calculation  of  the  shaded  histogram  was  about 
16°  («s  cos  *  0.961)  .  Considering  only  these  two  angles,  we  might 
therefore  expect  the  area  under  the  curve  to  be  greater  than  the  area 
of  the  shaded  histogram.  On  the  other  hand,  the  lead  collimator  scatters 
radiation  less  efficiently  than  the  steel  liner  of  the  collimator  in 
Figure  1.  It  thus  appears  that  the  larger  angle  of  acceptance  in  the 
Cormack  and  Johns  calculation  is  approximately  compensated  by  the 
relatively  poorer  scattering  characteristics  of  lead. 
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IV.  CONCLUSION 

The  data  presented  in  this  report  show  that  roughly  17%  of  the 
photons  which  emerge  with  their  initial  energies  from  the  collimator  in 
Figure  1  have  penetrated  some  portion  of  the  lip  of  the  collimator.  This 
number  could  probably  be  reduced  somewhat  if  the  steel  liner  of  the 
collimator  were  replaced  by  lead. 

It  was  also  observed  that  for  every  100  photons  which  are  emitted 
from  the  collimator  without  energy  loss,  approximately  32  photons  are 
scattered  by  the  collimator  wall  and  emerge  along  directions  which  cannot 
"see"  the  source  directly.  Most  of  these  photons  have  energies  in  excess 
of  0.5  MeV  with  the  most  probable  energy  being  about  1  MeV  .  This 
result  is  in  qualitative  agreement  with  earlier  calculations  published 

(5) 

by  Cormack  and  Johns  .  These  collimator- scattered  photons  were  not 

(2) 

observed  in  the  experimental  measurements  of  Costrell  mainly  because 
of  the  geometrical  arrangement  in  which  his  measurements  were  made. 
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Figure  1.  Source  in  collimator.  The  position  of  an  auxiliary  collimator 
which  was  placed  on  the  main  collimator  when  the  experimental  source  was 
in  use  is  indicated  by  the  dashed  lines  at  the  top  of  the  figure.  The 
auxiliary  collimator  was  not  included  in  the  computational  monel. 
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Figure  3.  Comparison  with  spectrum  calculated  by  Preiss  for  an  unco 
limated  source.  Error  bars  indicate  standard  deviation  of  our  Monte 
Carlo  results. 
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Spectrum  of  radiation  emerging  from  the  collimator  in  the 
interval  0.961  <,  u  <  0.996  . 


direction  interval  0  £  u  <  0.961 
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Figure  7.  Comparison  of  present  results  with  calculation  by  Cormack  and 
Johns  of  the  spectrum  of  radiation  singly  scattered  from  the  walls  of  a 
collimator. 
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APPENDIX  A 


FILTERED  SPECTRA 


For  use  in  experimental  work  at  the  University  of  Illinois,  it 
was  necessary  to  modify  the  collimator  emergent  spectrum  of  Figure  2 
by  assuming  that  radiation  from  the  collimator  subsequently  passed 
perpendicularly  through  a  thin  lead  sheet.  A  simple  filtered  spectrum 
was  calculated  in  which  scattering  by  the  lead  sheet  was  neglected. 

The  results  are  given  in  Table  Al.  In  the  column  for  zero  thickness, 
the  spectrum  was  normalized  to  one  photon  with  energy  greater  than 
0.1  MeV  .  In  the  other  columns,  for  a  given  energy  interval,  the  num¬ 
ber  to  the  left  of  the  slash  was  obtained  by  multiplying  the  corres¬ 
ponding  number  from  the  unfiltered  spectrum  by  an  exponential  factor. 
The  attenuation  coefficient  was  evaluated  at  the  mid-point  energy  of 
the  interval.  The  number  to  the  right  of  the  slash  resulted  from 
renormalizing  the  filtered  spectrum  to  one  photon  with  energy  greater 
than  0. 1  MeV  . 

The  spectra  in  Figures  2  and  5  were  added  together  and  filtered 
as  described  above  for  the  spectrum  in  Figure  2.  The  resulting 
spectra  for  this  case  are  presented  in  Table  A2. 
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TABLE  A1 


Number  of  photons  per  0.1  MeV  for  0.996  £  u  £1.0 
which  are  transmitted  by  lead  sheets  of  different 
thicknesses. 


Filter 
Thickness 
(inches) 

'V. 

Energy 
Interval 
(MeV) 

1.4  >  E  a  1.3 
1.3  >  E  £  1.2 
1.2  >  E  a  l.i 
l.i  >  E  a  l.o 
l.o  >  e  a  0.9 
0.9  >  e  a  o.s 
0.8  >  E  a  0.7 
0.7  >  E  a  0.6 
0.6  >  E  a  0.5 
0.5  >  r.  a  0.4 
0.4  >  E  a  0.3 
0.3  >  E  a  0.2 


0 

1/8 

1/4 

0.436 

0.358/0.474 

0.294/0.491 

0.003 

0.002/0.003 

0.002/0.003 

0.403 

0.324/0.429 

0.261/0.435 

0.011 

0.009/0.012 

0.007/0.012 

0.005 

0.004/0.005 

0.003/0.005 

0.010 

0.007/0.010 

0.005/0.009 

0.012 

0.009/0.012 

0.006/0.011 

0.011 

j  0.007/0.010 

0.005/0.008 

0.017 

j  0.010/0.014 

0.006/0.011 

0.020 

0.010/0.013 

0.005/0.009 

0.026 

1  0.010/0.013 

0.003/0.006 

0.032 

1  0.004/0.005 

0.001/0.001 
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TABLE  A 2 

Number  of  photons  per  0.1  MeV  for  0.961  £  uz  <1.0 
which  are  transmitted  by  lead  sheets  of  different 
thicknesses. 


Filter 

X.  Thickness 

(inches) 

Energy  'N. 

Interval 

(MeV) 

0 

1/8 

i 

1/4 

1.4  >  E  £  1.3 

0.390 

0.320/0.426 

0.263/0.443 

1.3  >  E  £  1.2 

0.016 

0.013/0.018 

0.011/0.018 

1.2  >  E  £  1.1 

0.379 

0.305/0.406 

0.245/0.413 

1. 1  >  E  £  1.0 

0.038 

0.030/0.040 

0.024/0.040 

1.0  >  E  £  0.9 

0.018 

0.014/0.018 

0.011/0.018 

0.9  >  E  £  0.8 

0.015 

0.011/0.015 

0.008/0.014 

0.8  >  E  £  0.7 

0.014 

0.010/0.013 

0.007/0.012 

0.7  >  E  £  0.6 

0.016 

0.011/0.015 

0.008/0.013 

0.6  >  E  £  0.5 

0.021 

0.013/0.017 

0.008/0.014 

0.5  >  E  £  0.4 

0.020 

0.010/0.014 

0.005/0.009 

0.4  >  E  £  0.3 

0.024 

0.009/0.012 

0.003/0.005 

0.3  >  E  £  0.2 

0.032 

0.004/C.  005 

0.001/0.001 

0.2  >  E  £  0.1 

0.016 

0 

0 

Total 

1.00 

0.750/1.00 

0.594/1.00 
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APPENDIX  B 

MASS  ATTENUATION  COEFFICIENTS 


The  mass  attenuation  coefficient  data  used  in  the  calculation  are 

listed  in  Table  Bl.  The  coefficients  for  iron  and  lead  are  from  tabu- 

171 

lations  of  Hubbell  and  Berger.  Since  Hubbell  and  Berger  did  not 
include  data  for  cobalt  in  their  tabulations,  data  for  this  element 
were  generated  by  a  separate  calculation  as  follows: 

The  total  Compton  interaction  cross  section  was  calculated  using 
a  tabulation  of  the  Klein-Nishina  formula.  The  cross  sections  for  three 
other  processes,  viz. ,  photoelectric  effect,  pair  production  in  the 
field  of  the  nucleus,  and  pair  production  in  the  field  of  an  electron, 
were  determined  by  an  interpolation  procedure  using  microscopic  cross 
section  data  for  iron  and  copper.  Three  priman  steps  were  involved  in 
the  interpolation. 

1.  A  dominant  functional  dependence  on  the  atomic  number 
/  was  removed  from  the  cross  section  data  for  iron 
and  copper  bv  division. 

2.  For  fixed  energy,  «»  ’meat  interpolation  was  performed 
on  the  data  for  iron  and  copper  to  get  corresponding 
data  for  cobalt. 

3.  The  cobalt  data  resulting  from  Step  2  were  multiplied 
by  the  function  of  Z  removed  in  Step  1. 

4.  5 

The  functional  dependences  assumed  were  Z  '  for  the  photoelectric 
2 

effect,  Z  for  pair  production  in  the  field  of  the  nucleus,  and  Z 
for  pair  production  in  the  field  of  an  electron. 


23 


TABLE  B1 


Mass  Attenuation  Coefficients 


(cm2/g)* 


Energy  (MeV) 

Iron 

Cobalt 

Lead 

0.010 

172.0 

187.0 

132.0 

0.015 

55.  7 

60.8 

112.0 

0.020 

25. 1 

27.4 

83.4 

0..  030 

7.87 

8.67 

27.9 

0.040 

3.46 

3.82 

13.0 

0.050 

1.84 

2.02 

7.17 

0.060 

1. 13 

1.24 

4.47 

0.080 

0.550 

0.596 

2.12 

0.100 

0.342 

0. 365 

5.62 

0.150 

0.184 

0.189 

1.99 

0.200 

0.139 

0.140 

0.969 

0.300 

0.107 

0. 106 

0.385 

0.400 

0.0921 

0.0911 

0.221 

0.500 

0.0829 

0.0818 

0.154 

0.600 

0.0761 

0.0751 

0.120 

0.800 

0.0664 

0.0655 

0.0856 

1.000 

0.0596 

0.0587 

0.0689 

1.500 

0. 0486 

0.0479 

0.0509 

2.000 

0.0425 

0.0419 

0.0450 

3.000 

0.0361 

0.0357 

0.0415 

4.000 

0.0331 

0.0328 

0.0415 

5.000 

0.0315 

0.0313 

0.0424 

6.000 

0.0305 

0.0304 

0.0434 

8.000 

0.0299 

0.0299 

0.0460 

10.000 

0. 0299 

0.0300 

0.0487 

3  3 

The  densities  assumed  were  37-86  g/cm  for  iron,  8.71  g/cm 
for  cobalt,  and  11.43  g/cm  for  lead.  The  density  for  lead 
should  have  been  11.34  g/cm3.  The  influence  of  this  mistake 
on  the  final  results  is  expected  to  be  negligible. 
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APPENDIX  C 

DESCRIPTION  OF  FORTRAN  PROGRAM  C0LLIMAT0R 

C, 1  General  Remarks 

C0LLIMAT0R  performs  a  Monte  Carlo  calculation  of  gamma  ray  transport 
in  a  finite,  multi-region  medium  having  cylindrical  symmetry.  Boundaries 
of  regions  are  defined  by  the  walls  of  concentric  cylinders  and  by  planes 
passing  through  these  cylinders  perpendicular  to  their  axis.  A  particu¬ 
lar  region  is  identified  by  a  radial  index  NR  and  a  vertical  index 
NZ  .  A  radial  index  NR  indicates  that  the  region  is  between  concentric 
cylinders  numbered  NR  and  NR  1  .  The  innermost  cylinder  is  assigned 

the  number  two.  Similarly,  a  vertical  index  NZ  means  that  the  region 
l.es  between  planes  numbered  NZ  and  NZ  +  1  . 

In  addition,  each  region  is  assigned  three  other  indexes:  MAT, 

VAT2,  ,.n<!  N'SPL.  MAT  '  <’ 1  ed  pseudo-nnterial  index.  If  it  is  posi¬ 

tive,  it  identifies  the  matenal  whose  attenuation  coefficient  is  to  be 
used  f«>r  s'mpling  path  lengths  ir  the  region  (sec  Section  II).  If  it 
is  zero,  the  region  -  ioid.  Tie  value  -1  indicates  that  the  region 
is  blr.  :>id  the  >  al  e  -2  lands  tl  re  ua>i  ,  v  detector. 

The  index  MAT2  i  either  positive  or  zero.  It  identifies  the 
actual  material  in  the  region.  If  for  a  particular  region  its  value 
is  different  from  the  value  of  MAT  for  the  region,  the  fictitious 
interaction  described  in  Section  II  is  introduced  as  an  alternative  to 
photoelectric  absorption  or  Compton  scattering.  MAT2  is  assigned  the 
value  zero  for  detector  regions,  black  regions,  and  voids. 

The  third  index,  NSPL,  is  a  splitting  parameter.  Splitting  was 
not  used  in  the  calculation  for  this  report,  but  the  program  was  written 
to  allow  for  this  possibility.  If  the  NSFL  value  for  a  region  is  zero, 
then  no  splitting  occurs  in  that  region.  If  the  value  of  NSPL  is  greater 
than  zero,  then  upon  interacting  in  the  region,  a  photon  is  split  into 
NSPL  new  photons.  The  statistical  weight  assigned  to  these  new  photons 
is  the  statistical  weight  of  the  original  photon  divided  by  NSPL.  When 
a  splitting  occurs,  the  radial  and  vertical  index  of  the  region,  the 
coordinates  of  the  interaction  point,  the  statistical  weight,  the  photon 
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energy,  and  the  direction  cosines  of  the  photon  are  stored.  The  history 
for  each  newly  generated  photon  is  then  followed  in  turn  until  the 
histories  of  all  the  photons  have  been  completed.  A  photon  which  results 
from  a  splitting  may  itself  be  split  if  it  has  a  subsequent  interaction 
in  a  region  where  NSPL  is  different  from  zero.  A  logic  diagram  for  a 
C0LLIMAT0R  history  with  splitting  included  is  shown  in  Figure  Cl.  In 
the  diagram,  ISPL  is  an  index  which  counts  the  number  o<  interaction 
points  where  splitting  has  occured. 

C. 2  Description  of  Individual  Subroutines 

C0LLIMAT0R  (Main  Program) 

Subroutine  called:  LEARN,  DATEX,  TRACK,  TEACH. 

The  program  consists  of  a  nest  of  five  loops.  The  index  NZS0 
specifies  the  z  coordinate  of  the  source  point,  NRS0  the  x  coordinate 
cl  the  source  point,  (the  y  coordinate  oi  the  source  point  is  always 
zero),  NES0  specified  the  source  energy,  and  NSPE  counts  the  number  of 
histories  completed  for  each  source  point,  source  energy  combination. 

The  product  NSEG  *■  NESX  x  NPSX  *  NZSX  *  LP  is  the  total  number  of 
histories  to  be  done. 

Subroutine  LEARN 

Subroutine  called:  none. 

All  input  data  are  read  by  this  program.  After  the  data  have  been 
read,  they  are  immediately  printed  under  headings  which  give  either  the 
variable  name  or  a  brief  description  of  the  variable. 

The  first  card  read  contains  72  alphanumeric  characters  and  when 
printed  at  the  beginning  of  the  calculation  serves  to  identify  the  run. 
Input  variables  are  listed  and  defined  b^low  in  the  order  in  which  they 
occur  in  the  subroutine. 

NRAN:  The  number  of  random  numbers  to  be  rejected  before  starting  the 
calculation. 

NEAX:  The  number  of  energy  intervals  to  be  used  in  the  calculation  of 
spectra. 

The  number  of  angular  intervals  to  be  used  in  calculating  the 
directional  distribution. 


NTHX: 
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MXAX:  The  number  of  energies  for  which  tabulated  cross  section  data  are 
to  be  read. 

MATX:  The  number  of  materials  m  the  collimator.  In  the  present  cal¬ 
culation,  there  were  three;  iron,  cobalt,  and  lead, 

NRAX:  The  number  of  radial  boundaries. 

NZAX:  The  number  of  horizontal  boundaries. 

NZTP:  The  index  of  the  horizontal  plane  bounding  the  top  of  the 

collimator. 

NRWL:  The  index  of  the  radial  boundary  of  the  air  column. 

NZBT:  The  index  of  the  horizontal  plane  bounding  the  bottom  of  the  air 
column. 

IHIX:  This  integer  must  be  divisible  by  the  product  NESX  *•'  NRSX  *  NZSX 

LP:  The  product  of  LP  and  IKIX  determines  the  number  of  histories  to 

be  done. 

NESX:  The  number  of  source  energies.  Each  source  energy  is  regarded  as 
equally  likely, 

NRSX:  The  number  of  x-coordinates  for  source  points. 

NZSX:  The  number  of  z-coordmates  for  source  points.  The  number  of 

source  points  is  NZSX  "  NRSX  . 

NR1:  The  radial  index  of  the  source  region. 

NZ1:  The  vertical  index  of  the  source  region, 

NASX:  This  index  will  be  discussed  later  with  the  variable  CTH0  and 
S0ANG  , 

MAT, 

MAT 2:  These  variables  were  discussed  m  the  general  remarks  at  the 

beginning  of  this  Appendix. 

ES0:  A  list  of  source  energies. 

ZH:  A  list  of  horizontal  boundaries. 

RH0:  A  list  of  cylindrical  boundaries.  The  first  number  in  this  list 

should  be  zero, 

DEN:  A  list  of  densities  for  the  materials  m  the  collimator, 

EB:  A  list  of  energies  for  which  cross  section  data  are  to  be  read. 

This  list  is  read  for  each  material  bu*  must  be  the  same  for 
each  material, 

XSECB:  This  is  a  list  of  interaction  coefficients.  If  the  indices  speci 
fied  for  this  variable  are  fM.l.NMAT),  then  we  have  the  Compton 
interaction  coefficient  for  the  M'th  energy,  and  the  material 


28 


with  index  NMAT.  If  the  index  1  is  replaced  by  2,  then  we  have  the 
total  attenuation  coefficient.  Note  that  if  the  densities  were 
read  in  the  order  iron,  cobalt,  and  lead,  then  the  interaction  co¬ 
efficients  must  be  read  in  the  same  order. 

DUM1, 

DUM2 , 

DUM3 , 

DUM4:  These  are  variables  which  are  used  to  skip  interaction  coefficient 
data  on  the  input  cards  which  are  not  used  in  the  calculation. 

CTH:  A  list  of  lower  limits  for  the  intervals  in  which  y  is  calssl- 

z 

fied  for  the  directional  distribution. 

ECLS:  A  list  of  lower  limits  for  the  energy  intervals  used  in  the  cal¬ 
culation  of  spectra. 

NSPL:  The  splitting  parameters  were  discussed  in  the  general  remarks 

at  the  beginning  of  this  appendix.  They  should  not  be  given  the 
value  one  since  this  accomplishes  the  same  thing  as  the  value 
zero  but  increases  the  amount  of  bookkeeping  required. 

CTH0, 

S0ANG , 

NASX:  The  program  assumes  that  each  source  point  emits  radiation  iso¬ 
tropically.  However,  it  may  be  desirable  to  sample  more  photons 
in  a  certain  range  of  directions  and  give  them  smaller  weight 
while  sampling  fewer  in  other  directions  and  giving  them  larger 
weight.  These  variables  define  a  cumulative  probability  dis¬ 
tribution  at  NASX  values  of  the  direction  cosine  relative  to  the 


z-axls.  The  probability  distribution  is  assumed  to  be  in  histo¬ 
gram  form.  To  sample  all  directions  with  equal  probability,  we 
read  in  the  following  table: 


The  values  of  S0ANG  must  be  in  increasing  order. 

IWCT:  For  each  material  in  the  collimator,  this  variable  lists  the 

index  of  one  of  the  energies  in  the  list  ECLS.  We  are  thus  able 
to  specify  different  cutoff  energies  in  different  materials  in 
the  collimator. 


CTH0  S0ANG 
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Subroutine  DATEX 

Subroutines  called:  RAM2B,  TABIN. 

The  tasks  of  this  subroutine  are  listed  below  in  the  order  in  which 
they  are  done. 

1.  Rejects  NRAN  random  numbers. 

2.  Re-arranges  the  tables  of  interaction  coefficients  so  that 
the  list  EB  is  in  descending  order  and  multiplies  the  inter¬ 
action  coefficients  for  each  material  by  the  density.  Since 
log-log  interpolation  is  used,  it  takes  logarithms  of  EB  and 
XSECB. 

3.  Prepares  an  expanded  table  of  interaction  coefficients  for 
300  energies  between  zero  and  1.5  MeV.  The  energy  interval 
0.005  MeV  is  used  over  the  entire  range.  Hereafter,  when  an 
interaction  coefficient  is  desired  for  an  energy  E  ,  rather 
than  interpolate,  we  compute  an  index  using  the  formula  NE  = 
200.0  *  E  +  G.5  and  use  the  entry  in  these  expanded  tables 
corresponding  to  this  index. 

4.  A  table  of  sines  and  cosines  is  generated  for  360  angles  start¬ 
ing  at  0.5°  and  proceeding  in  steps  of  1°  up  to  and  including 
359.5°  .  The  sine  and  cosine  for  a  random  angle  between  0° 
and  360°  are  then  selected  by  computing  an  index  according 

to  the  formula  IPH  =s  360.0  *  RAN  ,  where  RAN  is  a  random 
number,  and  then  choosing  the  sine  and  cosine  in  this  table 
corresponding  to  the  index  IPH  +  1  . 

5.  A  list  RH02  of  the  squares  of  the  radii  of  the  cylindrical 
boundaries  is  generated. 

6.  Lists  of  source  point  coordinates  RS0  and  ZS0  are  generated. 

7.  A  list  PART  of  floating  point  values  of  NSPL  is  generated. 

8.  Variables  in  which  scores  are  accumulated,  SPEC  and  SPEC2, 
are  initially  set  equal  to  zero. 

9.  The  cutoff  energy  for  each  material  is  expressed  in  terms  of 
its  Compton  wavelength. 

Subroutine  TRACK 

Subroutines  called:  RAM2B ,  CHECK,  VACU0,  C0MPT2,  GRADE. 


This  subroutine  performs  the  following  tasks: 
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1.  Source  position,  energy,  and  direction  are  selected  to  start 
a  photon  history.  The  energy  of  the  photon  is  expressed  in 
terms  of  the  Compton  wavelength. 

2.  Distances  between  interaction  points  are  sampled. 

3.  The  type  of  interaction  is  selected. 

4.  If  splitting  is  used,  the  necessary  bookkeeping  is  performed 
by  this  subroutine. 

Subroutine  CHECK 

Subroutines  called:  none. 

Without  going  into  very  much  detail,  we  may  outline  the  operation 
of  the  subroutine  by  the  following  steps: 

1.  It  checks  to  see  if  a  horizontal  boundary  has  been  crossed. 

If  it  has,  it  then  checks  to  see  if  a  cylindrical  boundary 
was  crossed  first. 

2.  If  a  horizontal  boundary  has  not  been  crossed,  it  checks  to 
see  if  a  cylindrical  boundary  has  been  crossed.  If  not,  it 
sets  the  index  NSCT  =  1  and  returns  to  TRACK. 

3.  If  a  boundary  is  crossed,  it  compares  the  index  MAT  for  the 
new  region  with  MAT  for  the  region  wnere  the  photon  started. 

If  MAT  is-  the  ->ame  for  both  regions,  then  NR  and  NZ,  the 
radial  and  vertical  indices  (see  the  general  remarks  at  the 
beginning  of  this  Appendix)  are  set  equal  to  values  appropriate 
tc  the  new  region.  It  then  checks  as  in  steps  one  and  two  to 
see  if  the  photon  actually  stops  in  the  new  region  or  crosses 
into  still  another  region.  When  a  photon  has  stopped  in  a 
region,  the  subroutine  sets  NSCT  =  1  and  returns  to  TRACK. 

4.  If  MAT  for  the  region  where  the  photon  started  and  MAT  for 
the  new  region  are  different,  the  coordinates  of  the  cross¬ 
ing  point  are  computed.  It  sets  NSCT  =  0  and  returns  to 
TRACK. 

5.  If  MAT  for  the  new  region  is  negative,  the  subroutine  does 
not  compute  crossing  point  coordinates  but  immediately  returns 
to  TRACK. 

Subroutine  VACU0 

Subroutines  called:  none. 


31 


Control  is  transferred  to  this  subroutine  when  the  photon  enters 
the  air  column.  The  coordinates  of  the  point  where  the  photon  crossed 
into  the  air  are  given.  The  subroutine  moves  the  photon  along  its 
current  trajectory  until  the  photon  leaves  the  collimator  or  re-enters 
some  part  of  the  collimator.  If  the  photon  re-enters  the  collimator 
or  source  capsule,  the  coordinates  of  the  re-entry  point  are  computed. 

Subroutine  C0MPT2 

Subroutine  called:  RAM2B. 

(g) 

This  subroutine  uses  the  method  of  Kahn  to  sampre  a  new  photon 
direction  and  energy  from  the  Klein-Nishina  distribution. 

Subroutine  GRADE 

Subroutines  called:  none. 

When  a  photon  enters  a  region  designated  as  a  detector,  the  sub¬ 
routine  adds  the  statistical  weight  of  the  photon  to  SPEC  for  the 
appropriate  energy,  angular  interval.  The  squares  of  individual  scores 
are  cumulated  in  SPEC2. 

Subroutine  TEACH 

Subroutine  called:  none. 

This  subroutine  normalizes  and  prints  the  data.  It  also  computes 
the  fractional  standard  deviation  of  all  results.  In  addition,  the 
spectrum  which  is  obtained  by  summing  over  all  direction  intervals  and 
the  directional  distribution  summed  over  all  energy  intervals  are 
printed. 

Liberal  use  is  made  of  printed  headings  over  all  output  data. 

This  makes  the  output  fairly  self-explanatory. 

Subroutine  TABIN 

Subroutines  called:  none. 

This  subroutine  performs  a  three-point  interpolation  simultaneously 
on  each  of  UMAX  functions  of  a  single  independent  variable.  XB  is  the 
list  of  independent  variable  values  for  which  the  functional  values  FB 
are  given,  and  must  be  stored  in  descending  order.  NMAX  is  the  number 
of  values  XB  in  the  list  and  must  be  i  3  .  FX  contains  the  list  of 
functional  values  interpolated  at  X  . 
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The  index  NTABIN  must  be  unity  the  first  time  the  subroutine  is 
called.  The  subroutine  then  computes  and  stores  certain  differences 
and  products  which  will  be  used  over  and  over.  Once  this  has  been  done, 
j.TABIN  is  set  equal  to  two  so  that  these  differences  need  not  be  com¬ 
puted  again  on  subsequent  calls. 

Functional  Subroutine  RAM2B 

This  function  is  used  to  generate  a  sequence  of  random  numbers 

uniformly  distributed  between  zero  and  one.  It  is  written  in  machine 

language.  The  multiplicative  congruential  method  is  used.  The  mul- 
15  35 

tiplier  is  5  and  the  modulus  is  2  .  The  first  number  in  the 

sequence  is  always  [5*^  (2*^  -  1)]  (mod  2^)  .  Thus,  to  change  the 
initial  random  number,  one  rejects  several  random  numbers  before  begin¬ 
ning  the  calculation. 

C. 3  Program  Listings  and  Sample  Output 

Listings  of  the  various  FORTRAN  programs  are  given  on  the  following 
pages.  The  computer  output  for  the  calculation  described  in  this  report 
is  reproduced  following  the  listings. 
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$  fastran 

C  COLLIMATOR  4-20-67 

C  CALCULATES  THE  ENERGY-ANGLE  SPECTRUM  OF  RADIATION  EMERGING 

C  FROM  A  COLLIMATED  SOURCE  OF  GAMMA  RADIATION 

DIMENSION  CPH ( 360 ) iCTh(20) *DEN(4) , EB(25) *ECLS(30) *MAT( 10.10)' 

1  MAT2( 10* 10)  *NEREM(50) *  NPREM ( 50 )  *  NRREM ( 50 ) ,NSPL(  10* 10)  » 

2  NZREM ( 50 )  , PARTI  10*10)  *RHO( 10)  ♦ RH02 ( 10) *RSO(  10)  ♦ SPEC (20 *30)  * 

3  SPEC2 ( 20  *  30  >  ,SPH(360)  *UC(300*4) ,UTI300.4) ,UXREM(50)  * 

4  UYREM ( 50 )  , UZREM (50)  , WREM (50)  . WTREM ( 50 )  *XREM(50)  *XSEC(2*4  )  * 

5  XSECBI 25*2 i 4 ) « YREMI 50 ) . ZH ( 1 0 ) *ZREM<50 ) *ZSO( 1 0 ) 

DIMENSION  ES0I2) *CThO(  10)  *SOANG(  10)  *  IWCTI4) *  WMAX I  4 ) 

COMMON  CPH*  CTH  *DEN , EB, ECLS , ESO ,  1 H I  X  *  MAT , MAT2 , MATX * MXAX * 

1  NE*  NEAX*NEREM,NES0*NESX*NPREM*NR*NR1  *  NRAN*  NRAX  *  NRREM  *  NRSO  * 

2  NRSX  *  NRWL  « NSCT  *  NSEG , NS PL  *NTHX«NZ*NZ1  ,  NZAX * NZBT * NZREM , NZSO * 

3  NZ  SX , NZTP , PART , RAN , RHO . RH02 . RSO  *  S , SPEC  *  SPEC2 *  SPH  *  SUZ  *  UC  *  UT  * 

4  UX  *  UXREM , UY , UYREM , UZ  *  UZREM , W ♦ WMAX  *  WREM , WT ♦ WTREM , X  *  XN , XREM  , 

5  XSEC . XSECB  *  Y  *  YN . YREM . Z  *  ZH  *  ZN  *  ZREM , ZSO 
COMMON  CTHO*  IWCT  * NASX * SOANG * LP 

CALL  LEARN 
CALL  DATEX 
L  =  0 

5  L=L+ 1 
NZSO=  0 

10  NZSO=  NZSO+ 1 
NRSO^O 

20  NRSO=  NRSO+ 1 
NESO= 0 

30  NESO=  NESO+ 1 
NSPE=0 

40  NSPE=NSPE+1 
CALL  TRACK 

I F ( NS PE— NS EG ) 40  *  50  *  50 
50  IF CNESO-NESX J30.60 *60 
60  I  F  (  NRSO-*NRSX  >20*70*70 
70  IF(NZSO-NZSX) 10*80*80 
80  IFIL-LP)5. 90*90 
90  CALL  TEACH 
CALL  SYSTEM 
END 


$  FASTRAN 

C  SUBROUTINE  LEARN  4-20-67 

C  READS  ALL  INPUT  DA  i‘A 

SUBROUTINE  LEARN 

DIMENSION  CPH ( 360 )  »CTh( 20 )  *  DENI  4 )  ,EB (25 )  * ECLS ( 30  )  *  MAT  I  10*10)* 

1  MAT2 ( 10* 10 ) *NEREM( 50 ) * NPREM { 50 ) , NRREM ( 50 ) *NSPL (10*10)* 

2  NZREM (50) *PART(  10*10) *RhO( 10) «Rh02< 10) *RSO<  10)  *  SPEC (20.30)  * 

3  SPEC 2 (20, 30) , SPH (360) ,UC( 300,4 ) ,UT (300,4' , UXREM (50) , 

4  UYREM ( 50 ) , UZREM ( 50) * WREM <  50 )  * WTREM (50 )  ,XREM( 50 ) *XSEC(2*4 ) ♦ 

5  XSECB (25 *2* 4) , YREM (50)  , ZH ( 1 0 ) * ZREM (50)  *ZSO(  10) 

DIMENSION  ESO ( 2 ) *CThO( 10) , SOANG ( 10) * IWCT (4) .WMAX (4) 

COMMON  CPH.CTHiDEN, EB, ECLS, ESO *  I H I  X  *  MAT , MAT2 * MATX * MXAX * 

1  NE  *  NEAX , NEREM*  NESO , NESX  *  NPREM , NR  «  NR  1 , NRAN , NR AX  * NRREM , NRSO , 

2  NRSX*NRWL,NSCT*NSEG*NSPL,NTHX*NZ*NZ1 ,  NZAX.NZBT, NZREM, NZSO* 

3  NZ SX.NZTP, PART, RAN* RHO, RH02, RSO, S* SPEC, SPEC2, SPH, SUZ *UC*UT, 
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4  UX . UXREM i UY , UYREM , UZ i UZRCM , W , WMAX , WREM  *  WT , WTREM , X  *  XN . XREM , 

5  XSEC  «  XSECB  <Y,YNi YREM , Z i ZH  .  ZN , ZREM . ZSO 
COMMON  CThO  (.  IWCT.NASX.SOANG.LP 

W0T6*  10 
10  FORMAT ( 1 HO ) 

R 1T7. 20 
WOT& .  20 
20  FORMAT (72H 

1  ) 

W0T6 » 10 
W0T6.  30 

30  FORMAT (72H  NRAN  NEAX  NTHX  MXAX  MATX  NRAx  NZAX  NZTP  nrwl 

1NZST  I H 1 X  LP) 

R I T7 . 45 . NRAN , NEAX  .  NTHX  .  MX  AX  ,  MATX  . NRAX  .  NZAX . NZTP , NRwL , NZBT  ,  1 H 1 X , LP 
WQT6 «  45 .  NRAN  ,  NEAX  .NT  HX i MX AX  .  MATX  .  NRAX ♦  NZAX .  NZTP  » NRwL ,  NZBT i  I H I  X ♦ LP 
40  FORMAT! IB  16) 

45  FORMAT! 1016. 18* 14) 

W0T6  «  10 
W0T6  *  50 

50  FORMAT (36H  NESX  NRSX  NZSX  NR 1  NZ1  NASX ) 

R I T7 « 40  .  NESX  .  NRSX . NZSX  . NR  1 .NZl  .NASX 
W  0  T  6  « 40, NESX. NRSX .NZSX. NR  1  ,NZl  .NASX 
W0T6. 10 
W0T6 . 60 

60  FORMAT ! 65H  INDICES  TO  IDENTIFY  PSEUDO-MATERIAL  IN  EACH  REGION  OF  C 
1 OLL IMATOR) 

DO  70  NZ=1 .NZAX 

R IT7. 40. IMAT1NR.NZ) .NR=1 .NRAX) 

WOT  6 . 40 ,  !MAT!NR,NZ)  ,NR=l  ,NRAX) 

70  CONTINUE 
WO To .  1  o 
W0T6 . 80 

80  FORMAT ( 58H  INDICES  TO  IDENTIFY  MATERIAL  IN  EACH  REGION  OF  COLL  I  MAT 
1  OR) 

DO  90  NZ= 1 .NZAX 

RIT7.40, ! MAT2 ( NR . NZ ) .NR=1 .NRAX) 

WOT6  »  40.  ! MAT2 ( NR. NZ ) . NR= 1  , NRAX ) 

90  CONTINUE 
W0T6« 1 0 
WOT6. 120 

120  FORMAT! 16H  SOURCE  ENERGIES! 

RIT7. 130. ( ESQ ! NESO ) ,NESO= 1 ,NESX) 

W0T6. 130.  ! ESO ! NESO )  .NES0=1  ,NESX) 

130  FORMAT! 10F7.3) 

W0T6. 1 0 
WO  T  6 .  140 

140  FORMAT (37H  LOWER  BOUNDARIES  OF  VERTICAL  REGIONS) 

R1T7. 150. (ZH(NZ) .NZ=1 .NZAX) 

W0T6. 150. (ZH(NZ) ,NZ=1 .NZAX) 

150  FORMAT ( 9F8  »  2  ) 

WGT6 ♦  10 
W0T6. 160 

160  FORMAT  (35H  INNER  BOUNDARIES  OF  RADIAL  REGIONS) 

R I T  7 ,  150.  ! RhO ( NR )  ,NR=i  , NRAX ) 

W0T6. 150. (RHOINR) ,NR=1 , NRAX ) 

WO  T6 .  1  0 
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WGT6 4 170 

170  FORMAT ( 1 7H  DENSITIES^  G/CMJ) 

R I T7  »  180i  (DEN(NMAT)  ,NMAT=1  4 MATX) 

W0T6 i 180i (DEN(NMAT) ,NMAT=1 tMATX) 

180  FORMAT <  1P6E1 1 .3) 

W0T6  »  10 

DO  230  NMAT  =  1  ♦ MATX 
R1T7« 190 
W0T6  4  190 
190  FORMAT (72H 

1  ) 

W0T6f  200 

200  FORMAT  ( 24H0E’iERGY  COMPTON  TOTAL) 

DO  220  M= 1 iMXAX 

RIT7i210« EB { M )  «DUM1  ,XSECB(M,1  ,NMAT)  *  DUM2 i DUM3  *  DUM4 , XSECB (Mi2i NMAT ) 
V/0T6  »  2  1  0  t  EB  (  M  )  *  (  XSECB  (  M  »  I  «  NMAT  )  «  1  =  1  *  2  > 

210  FORMAT ( F7«  3  *  1 P7E 10*3) 

220  CONTINUE 
WO  T6  ♦  1  0 
230  CONTINUE 
W0T6 ♦ 240 

240  FORMAT ( 57H  LOWER  LIMITS  ON  COSINE  THETA  FOR  ANGULAR  CLASSIFICATION 
IS) 

R I T7 i 250 i (CTH(NTH) *NTh=1 »NTHX) 

W0T6 ♦ 250  « <CTH(NTH)  iNTh=1  *NTHX) 

250  FORMAT ( 8F9»5 ) 

W0T6«  10 
WO  T6 ♦ 260 

260  FORMAT (40H  LOWER  LIMITS  FOR  ENERGY  CLASSIFICATIONS) 

R I T7 i  1 30 i ( ECLS  <  NEC )  cNEC  =  l  ,N£AX> 

W0T6 i 130* (ECLS (NEC) ,NEC=1 ,NEAX) 

WO  T6  «  1  0 
W0T6  4  270 

27C  F  ORM  AT ( 2 1 H  SPLITTING  PARAMETERS) 

DO  280  NZ= 1 4  NZAX 

R I T7i 40  »  (NSPL ( NR i NZ )  i NR= 1  >  NR AX ) 

W0T6«40, (NSPLINR^NZ) .NR=1 ,NRAX) 

280  CONTINUE 
W0T6  ♦  10 
WO  T 6 4  290 

290  FORMAT ( 37H  SOURCE  SAMPLING  ANGULAR  DISTRIBUTION) 

R I T7  4  250  4  (CTHO(NASO)  iNASO=l  4NASX ) 

WOT64  250  4  ( CTHO ( NASO )  4NAS0=1  4  NASX ) 

R I T7  4  250  4  (SOANG(NASO)  4NAS0=1 4  NASX ) 

WOT  6  *  250  4  (SOANG(NASO) 4NAS0=1 4  NASX ) 

WOT64 10 
W0T6  v  300 

300  F  ORM  AT ( 22H  ENERGY  CUTOFF  INDICES) 

RIT744O4  ( IWCT(NMAT)  , NMAT =1  4  MATX ) 

WOT644O4  (  IWCT(NMA1  )  »NMAT=l  4  MATX ) 

W0T6< 10 

RETURN 

END 
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C  SUBROUTINE  DATEX 

C  EXPANDS  INPUT  DATA  FOR 

c  work  needed  for  the  ma 

SUBROUTINE  datex 
DIMENSION  CPH ( 360 ) <CTh 

1  MAT2( 10. 10) .NEREM(50 

2  NZREM  <  50 )  , PART (10*10 

3  SPEC2 ( 20  , 30 )  • SPH < 360 

4  UYREM ( 50 )  , UZREM ( 50 ) , WREMI50)  ,WTREM<50>  , XREM ( 50 ) .XSEC <2, 4  )  , 

5  XSECB ( 25 *2*4 ) i YREM( 50 )  ,  ZH ( 1 0 ) ♦ ZREM ( 50 )  ,ZSO(  I  0 ) 

DIMENSION  ES0(2>  *CThO( 10>  *  SOANG  <  10)  .  IWCT(4) , WMAX(4) 

COMMON  C»H  «  CTh  *DEN  ,  EB,  ECLS ,  ESO  ,  IHIXi MAT  «  MAT2  «  MATX  * MXAX , 

1  NE  ,  NEAX ,  NEREM  ,  N£SO , NESX i NPREM , NR , NR  1  , NR AN , NR AX , NRR£M , NRSO , 

2  NRSX i NRwL , NSCT  «  NSEG,  NS PL , NTHX , NZ , NZ 1  ,  NZAX  ,  NZBT , NZREM , NZSO , 

3  NZ  S  X  4  NZTP , PART , RAN , RhO  ,  RH02 , RSO , S , SPEC  *  SPEC2 . SPH  ,  SUZ , UC . UT , 
a  jx , wXREM , UY , UYREM  ,  UZ  .  UZREM , W , WMAX , WREM  ,  WT ,  WTREM  , X  »  XN , XREM  , 

5  XSEC , XSECB .Y.YN, YREM , Z  »  ZH , ZN  . ZREM , ZSO 

COMMON  CThO. 1WCT,NASX,S0ANG 
NT AS  I Nt l 

DO  10  NRA=1*NRAN 
10  RAN=RAM2B(0) 

MX AX 2= MX AX/2 
MM-MX  AX+ 1 
DO  15  M= 1 , MXAX2 
MM a MM “ 1 
DUMaEB ( M ) 

EB (M) =EB  5  MM ) 

E8(MM)=DUM 
DO  15  1=1*2 

DO  15  NMAT= 1 .MATX 
DUMr XSECB ( M  ,  1 ,  NMAT ) 

XSECB (M,  I  * NMAT ) =XSECB ( MM  «  I  « NMAT ) 

XSECB (MM, I , NMAT ) =DUM 
15  CONTINUE 

DO  30  M= 1 , MXAX 
DO  20  1=1,2 

DO  20  NMAT  = 1  , MATX 

XSECB(M, I ,NMAT)=XSECB(M, I , NMAT )*DEN( NMAT) 

20  XSECB (M, I , NMAT )=ELOG< XSECB <M, I ,NMAT> ) 

30  EB ( M ) =ELOG ( EB ( M ) ) 

MATX2=2*MATX 
DO  50  NE= 1 .300 
E  =  NE 

E=0.005*E 
E»ELOG( E ) 

CALL  T AB IN( NT AB I N, XSECB, EB, MXAX, MATX2  »E , XSEC) 

DO  40  NMAT  = 1  , MATX 

UC ( NE . NMAT ) =EXP<  XSEC (  1  , NMAT )  ) 

40  UT(NE.NMAT)=LXP(XSEC(2,NMAT) ) 

50  CONTINUE 
PH  =  -0  »5 

DO  60  I PH= 1 , 360 
PH  =  PH+ 1 «  0 
PHR=PH*0401 7453293 
CPH (  I  PH ) aCOS(PHR) 

60  SPH ( I DH ) =S I N ( PHR ) 


3-2-67 

TABLE  LOOK-UP  AND  DOES  OTHER  PRELIMINARY 
IN  PART  OF  THE  CALCULATION 

(20) , DEN ( 4 ) ,EB(25) , ECLS (30) , MAT ( 10,10). 

)  ,  NPREM  (  50  )  ,  NRREM  (  50  )  ,  NSPi-  (10,10), 

)  , RHO  < 10) i RH02 (  10) , RSO ( 10) , SPEC (20,30) , 

)  ,UC (300,4)  , UT ( 300 , 4 ) ,UXREM(50)  , 


O  6> 
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DO  70  NR  = 1  ,  NRAX 
70  RH02( NR) =RH0(NR)**2 
NSPTa  NRSX*NZSX 
NSPGa I  HI X/NSPT 
IF ( NSPG*NSPT- I  HI  X ) 80  * i 1 0 i 80 
80  W0T6 • 90 

90  FORMAT (62H0  NUMBER  OF  HISTORIES  NOT  DIVISIBLE  BY  NUMBER  OF  SOURCE 
i POINTS) 

100  CALL  SYSTEM 
110  FNRSX  =  NRSX 

DO  120  NR30=1  4  NRSX 
FNRSOcNRSO 

120  RSO(NRSO) =0.5* (SQRT<  {  ( FNRSX-FNRSO ) *Rh02 ( NR  1  ) +FNRS0#Rh02 ( NR 1  +  1  )  ) 

1  /FNRSX ) +SQRT < (  ( FNRSX-FNRSO+1 .0 ) *Rh02 ( NR  1  )  +  (FNRSO-l .0 )*Rh02(NRi 

2  +1) ) /FNRSX) ) 

FNZSX=NZSX 

DO  130  NZSO  = 1  4  NZSX 
FNZSO=NZSO 

130  ZSO ( NZSO ) = ( (2.0*FNZS0-1 .0)*ZH(NZ1+1 ) + ( 2 . 0* ( FNZSX-FNZSO ) + 1 .0) 

1  *ZH ( NZ 1  )  ) / <  2 • 0*FNZSX ) 

NSEG=  nspg/nesx 

I F ( NSEG*NESX-NSPG ) 140,160*140 
140  W0T6  *  1 50 

150  FORMAT < 86H0  NUMBER  OF  HISTORIES  FOR  EACH  SOURCE  POINT  NOT  D1VISIBL 
IE  BY  NUMBER  OF  SOURCE  ENERGIES) 

GO  TO  100 

160  DO  170  NR=1  *  NRAX 
DO  170  NZ= 1 ,NZAX 
170  PARTI NR,NZ) =NSPL(NR,NZ) 

DO  1 80  NTH= 1  *  NTHX 
DO  180  NEC= 1 , NEAX 
SPEC  (  NTH  *  NEC  )  <=0  »  0 
180  SPEC2(NTH,NEC) =0.0 
DO  190  NMATrl *  MATX 
NEC= I WCT ( NMAT ) 

190  WMAX  ( NMAT ) =0 • 5 1 1 /ECLS  <  NEC ) 

RETURN 

END 


FASTRAN 

SUBROUTINE  TRACK  3-13-67 

C  GENERATES  PHOTON  HISTORIES 

S". ROUTINE  TRACK 

DIMENSION  CPH ( 360 )  *CTh(20)  ,DEN(4)  ,EB<25)  , ECLS (30)  *  MAT  (  10,10), 

1  MAT  2(10,10), NEREM ( 50 )  , NPREM ( 50 )  , NRREM ( 50 ) , NSPL (10,10), 

2  NZREM ( 50 ) , PART (10,10), RhO (10), RH02 (10), RSO (10), SPEC (20,30), 

3  SPEC2 ( 20 , 30 ) ,SPH(360) , UC < 300,4) ,UT( 300,4) ,UXREM(50) , 

4  UYREM ( 50 )  , UZREM ( 50 )  , WREM ( 50 )  , WTREM ( 50 )  » XREM ( 50 )  *  XSEC (2,4), 

5  XSE CB ( 25 ,2,4) ,YREM(50)  , ZH ( 1 0 ) ,ZREM(50)  ,ZSO( 10) 

DIMENSION  ESO ( 2 )  ,CThO(  10)  ,SOANG(  10)  ,  IWCT(4)  ,WMAX(4) 

COMMON  CPH, CTh, DEN, EB, ECLS, ESO,  I H I  X , MAT , MAT2 , MATX , MXAX , 

1  ne» neax, NEREM, NESO.NESX, NPREM, NR, NR  1  , nran, NRAX, nrrem, NRSO, 

2  NRSX, NRWL,NSCT,NSEG, NSPL, NTHX, NZ, NZ1 ,  NZAX , NZBT , NZREM , NZSO , 

3  NZSX, NZTP, PART, RAN, RH0,RH02, RSO, S, SPEC, SPEC2,SPH,SUZ,UC,UT , 

4  UX , UXREM , UY , UYREM , UZ , UZREM , W , WMAX , WREM , WT , WTREM , X , XN , XREM , 
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5  X,S£  C i XSECB *  Y  *YNt YREM  *  Z  *  ZH  <  ZN  *  ZREM  «  ZSO 
COMMON  CThO* lWCT*NASX,SOANG 
DIMENSION  SUZREM(50) 

I SPL  =  0 
XapSO(NRSO) 

Y=0  .0 

Z«ZSO(NZSO) 

W=  0  «  5 1 1 /ESO ( NESO  > 

RANcRAM2B(C) 

DO  2  NAS0  =  2  « NASX 
NAaNASO 

IF(RAN-SOANG<NASO) >4*2*2 
2  CONTINUE 
4  RAN=R AM2B 1 0 ) 

DC  THs  CTHO  <  NA ) — CTHO ( NA— 1  ) 

UZ  =  CTHO  (  NA—  1  )+RAN*DCTH 

WT  =  DCTh/(  <S0ANG(NA)-S0ANG<NA-1  )  )*2.0  ) 
RAN=RAM2B ( 0 ) 

IPH=360.0*RAN 
SUZ»  SORT ( 1 . 0-UZ**2 ) 

UX=SUZ*CPH( IPH+i > 

UY=3UZ*SPH(  IPH+1  ) 

NRcNR  1 
NZaNZI 

NE  =  1  0  2  •  2/W+O  •  5 
10  NMAT=MAT(NR,NZ) 

RANsR AM2B  <  0 ) 

NEaNE 

Sn -ELOG ( RAN ) /UT ( NE , NMAT ) 

XN=UX*S+X 

YN=UY*S+Y 

ZN=UZ*S+Z 

CALL  CHECK 

X=XN 

Y=  YN 

Z*ZN 

NR=  NR 

NZ  =  I  ’ 

IF(Mm  .NR*NZ) >80*20*30 
20  CALL  VACUO 
NR  =  NR 
NZ  =  NZ 

IF <MAT( NR, NZ) >80*30,30 
30  IF(NSCT> 10* 10*40 
40  IF(NSPL(NR*NZ) >60*60*50 
50  I SPL= I SPL+ 1 

I F ( ISPL-50J55, 55*60 
55  WT  =  W  T /PART ( NR  *  NZ ) 

NPREM  5  I SPL  >  aNSPL ( NR , NZ ) -  1 

XREM ( I SPL ) =X 

YREM  <  I SPL ) a Y 

ZREM (  I SPL ) c  Z 

UXREM  <  ISPL)nUX 

UYREM  <  I SPL) eUY 

UZREM  <  I SPL ) tUZ 

SUZREM  <  I SPL) =SUZ 
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WREM  <  ISPL ) =W 
NEREM ( ISPL ) =NE 
WT REM (  I SPL  >  =WT 
NRREM  (  I  SPL  )  siNR 
NZREM ( I SPL ) sNZ 
60  NE=NE 

NMATaMAT (NR i NZ ) 

NM ATR=  MAT2 ( NR  ,  NZ ) 

0  I  FA  T«  UT ( NE , NMAT ) -UT ( NE  .  NMATR ) 
DENOMnD IFAT+UC ( NE  *  NMATR ) 
RAN=RAM28<0 ) 

I F  ( RAN-DENOM/UT ( NE i NMAT )  ) 65  ,  65  , 1 00 
65  RAN=RAM28(0) 

IF ( RAN— D IFAT/DENOM ) 1 0 , 70 , 70 
70  CALL  C0MPT2 

NMAT=MAT2(NR,NZ) 

I F ( Wu WMAX ( NMAT ) ) 10* 1  00  , 100 
80  IF <MAT(NR,NZ)+1 )90» 1 00 » 1 00 
90  CALL  GRADE 
100  I F ( ISPL) 140. 140* 1 10 
110  I F ( NPREM ( I SPL ) ) 120t 120'130 
120  I SPL= 1 SPL- 1 
GO  TO  100 

130  NPREM ( ISPL) =NPREM< ISPL1-1 
X=XREM( ISPL) 

Y  =  YRE  M  <  I SPL) 

Z  n ZRE  M (  ISPL) 

UXeUXREMf I SPL ) 

UYaUYREM ( ISPL) 

UZ^UZREM ( ISPL) 

SUZeSUZREMC ISPL) 

W*  WRE  M (  ISPL) 

NE=NEREM ( ISPL ) 

WT  a WTREM  < ISPL) 

NR^NRREM ( ISPL ) 

NZ  =  NZ  REM ( ISPL) 

GO  TO  60 
1  AO  RETURN 
END 


$  FASTRAN 

C  SUBROUTINE  CHECK  3-2-67 

C  CHECKS  TO  SEE  IF  BOUNDARY  HAS  BEEN  CROSSED 

SUBROUTINE  CHECK 

DIMENSION  CPH ( 360 ) »CTH(20) .DEN (4) ,EB(25) .ECLSC30) . MAT ( 10.10). 

1  MAT2I10. 10) .NEREM (50) .NPREM (50) .NRREM (50) . NSPL ( 1 0 ♦ 1 0 ) « 

2  NZREM (50) . PART ( 10.10) . RHO ( 10) .RH02< 10) *RSO( 10) . SPEC ( 20 . 30 ) . 

3  SPEC2 ( 20 . 30 ) ,SPH(360) .UC(300.4) ,UT (300.4) ,UXREM(50) . 

4  UYREM ( 50 ) . UZREM ( 50 ) , WREM (50) .WTREM (50) .XREM(50) .XSEC(2.4 ) . 

5  XSECB( 25.2.4 ) . YREM ( 50 ) .ZH( 10) .ZREM(50) .ZSO( 10) 

D I  MENS  ION  ESO( 2)  . CThO(  1 0 )  . SOANG (10)  ♦  I WCT ( 4 ) . WM AX ( 4  ) 

COMMON  CPH.CTH.DEN.EB.ECLS.ESO,  I H 1 X  »  MAT . MAT2 ♦ MATX  «  MX AX ♦ 

1  NE. NEAX.NEREM.NESO.NESX. NPREM, NR, NR1  , NR AN , NR AX , NRREM  «  NRSO , 

2  NRSX , NRWL  «  NSCT  «  NSEG, NSPL .NTHX.NZ.NZ1  ,  NZAX.NZBT , NZREM, NZSO, 

3  NZSX , NZTP , PART , RAN , RHO , RH02 , RSO , S , SPEC  »  SPEC2  ,  SPH , SUZ , UC . UT  c 
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4  UX , UXREM i UY  .  UYREM , UZ . UZREM , W . WMAX , U/REM . WT , WT REM  . X . XN , XREM , 

5  XS£C  .  XSECB  *  Y . YN . YREM . Z . ZH  *  2N .  ZREM , ZSO 
COMMON  CTHO. 1WCT* NASXiSOANG 

DR AD=  X*UX+Y#UY 
IF ( DRAD) 20C .10.10 
10  IF(UZ) 130. 100,20 
20  IF  ( ZN— ZH ( NZ+ 1  )  )  l 00 . 30  «  30 
30  SB=(ZH(NZ+1)-Z)/UZ 
XI =S8*UX+X 
Y l =S3*UY+Y 

40  IF(X1**2+Y1**2~RH02(NR+1 > >90.70*50 
50  NR=NR+l 

IF (MAT(NR.NZ)-MAT(NR-l  ,NZ)  ) 60  *  40 . 60 
60  I F ( MAT ( NR .  NZ ) ) 480  t65i 65 
65  SUZ2=UX**2+UY**2 

SSs ( SORT (ORAD**2-SUZ2*<  X**2  +  Y**2-PH02 (NR)  )  ) —DRAD ) /SUZ2 
ZN-SB#UZ+Z 
XN»SB*UX+X 
YN*SB*UY+Y 
NSCT  b O 
GO  TO  480 
70  NR=NR+ 1 
NZ=NZ+ 1 

IF(MAT(NRiN7)-MAT(NR-1 .  NZ  —  1 ) >80.20.80 
80  I F ( MAT ( NR »  NZ )  ) 480 . 85  .  85 
85  ZNrZH(NZ) 

XN=X1 
YN  =  Y  1 
NSCT -  0 
GO  TO  480 
90  NZ  =  NZ+ 1 

IF !MAT(NR.NZ )-MAT(NR.NZ-l ) ) 80 » 20 *80 
100  IF(XN**2+YN**2-RH02<NR+i  )  )  1  20 .  1 1 0 . 1  1 0 

110  NR=NR+ 1 

IF (MAT(NR,NZ )-MAT(NR-l ,NZ) >60.100*60 
120  . sw  CT  =  1 

CO  TO  480 

130  IFCZH  <NZ  )*-ZN>  100. 140. 140 
140  SB=CZH(NZ)-Z)/UZ 
X  1  =SB“'OX+X 
Y 1 =SB*UY+Y 

150  IF(X1**2+Y1**2-Rh02(NR+1  )  )  1 90  *  1 70  *  1 60 

160  NR=NR+ 1 

IF (MAT(NR.NZ)-MAT(NR-l , NZ) >60.150*60 
170  NR=NR+ 1 
NZ*NZ~1 

IF (MAT(NR.NZ)-MAT(NR-I , NZ+1 ) ) 180. 130* 180 
180  IF(MAT(NR.NZ) )480. 185. 185 
185  ZN=  ZH ( NZ+ 1  ) 

XN=  X  1 
YN  =  Y  1 
NSCT  x 0 
GO  TO  480 
190  NZ=NZ“ 1 

IF<MAT(NR.NZ)-MAT(NR,NZ+1  ))  180. 130. 180 
200  SUZ2=UX**2+UY**2 


41 


SM= -DR AD/SUZ2 
DS “S— SM 

IF (DS)220.220,210 
210  XM=UX*SM+X 
YMaUY *SM  +  Y 
ZM=UZ*SM+Z 
GO  TO  230 
220  XM=XN 
YMaYN 
ZM  =  ZN 

230  IF ( UZ ) 390  .  340  «  240 
240  IF  (ZM-ZH(i'lZ+l  )  >340.250.250 
250  SB=(ZH(NZ+1 )-Z)/UZ 
X 1 =SB*UX+X 
Y1 =SB*UY+Y 

260  IF ( XI **2+Yl **2-RH02(NR ) ) 320.290 . 270 
270  NZ=NZ+1 

IF(MAT(NR,NZ)-MAT(NR,NZ-1 ) >280.240.280 
280  IF(MAT(NR.NZ) >480.285.285 
285  ZN=ZH(NZ) 

XN=  X  1 
YN  =  Y1 
NSCT s  0 
GO  TO  480 
290  NR=NR— 1 

I F ( NR ) 300 . 300 .310 
300  NR=NR+1 

GO  TO  270 
310  NZ=NZ+1 

IF(MAT(NR.NZ)-MAT(NR+l  ,NZ-1  >  >280.240.280 
320  NR=NR— 1 

I F (MAT ( NR , NZ ) -MAT ( NR+ 1  , NZ  >  >  330 . 260 . 330 
330  IF (MAT (NR.NZ ) >480 .335.335 

335  SB=< -DR AD-SORT < DRAD**2— SUZ2* ( X**2+Y**2-RH02 < NR+ 1 ) ) ) )/SUZ2 
XN=SB*UX+X 
YN=SB*UY+Y 
ZN=SB*UZ+Z 
NSCT  s 0 
GO  TO  480 

340  IF ( XM**2+YM**2— RH02 (NR ) >380.350.350 
350  IF (DS) 370. 370,360 
360  X=XM 
Y=  YM 
Z-  ZM 

DR  AD a  X*UX+ Y#UY 
GO  TO  10 
370  XN=XM 
YNaYM 
ZN=ZM 
NSCT  c  1 
GO  TO  480 
380  NR=NR— 1 

IF(MAT(NR.NZ)-MAT(NR+1 ,NZ) >330.340.330 
390  IF ( ZH ( NZ 1—ZMJ340. 400.400 
400  SB=( ZH(NZ)-Z)/UZ 
XI =SB*UX+X 
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Y1  =S£3*UY+Y 

410  IF (XI **2+Yl**2-RH02(NR) >470.440,420 
420  NZ=NZ— 1 

I F  <  MAT ( NR  ,  NZ ) —MAT ( NR , NZ  + i  )  ) 430 , 390*430 
430  IF(MAT(NR.NZ) >480,435,435 
435  ZN=ZH(NZ+1 ) 

XN  =  X1 
YN=Y1 
NSCTcO 
CO  TO  480 
440  NR=NR— 1 

IF (NR) 450* 450,460 
450  NR=NR+ 1 

GO  TO  420 
460  NZ=NZ— 1 

1F(MAT(NR,NZ)-MAT(NR+1 ,  NZ+1 ) >430.390,430 
470  NR=NR-1 

IF (MAT ( NR  ,  NZ ) ^MAT ( NR+ 1 ,NZ)  >330.410*330 
480  RETURN 
END 


$  FASTRAN 

C  SUBROUTINE  VACUO  3-2-67 

C  MOVES  PHOTON  FROM  POSITION  IN  HOLE  TO  WALL  OR  FLOOR  OF  COLLI 

C  OR  TO  DE' ECTOR 

SUBROUTINE  VACUO 

DIMENSION  CPH( 360 )  .  CTh<20 )  *DEN(4  >  ,EB(25)  ,ECLS(30  >  .MAT (10.10) 

1  MAT  2 (  10.10)  .NEREM(SO)  ♦ NPREM ( 50 )  , NRREM ( 50 ) ,NSPL( 10.10). 

2  NZREM ( 50 ) . PART (10.10). RhO (10). RH02 (10). RSO (10). SPEC (20.30) 

3  SPEC2 ( 20  «  30 )  ,SPh(360) ,UC(300.4) .UT(300.4) .UXREM(50) , 

4  JYREM ( 50 )  . UZREM ( 50 ) , WREM ( 50 )  . WTREM ( 50 )  . XREM ( 50 ) *  XSEC (2.4). 

5  XSECB ( 25  *2*4)  . YREM ( 50 )  . ZH ( 1 0 )  .ZREM(50)  .ZSO(  10) 

DIMENSION  ESO(2)  .CThO(  10)  .SOANG( 10)  .  IWCT(4) .WMAX(4) 

COMMON  CPH.CTH.DEN.EB.ECLS.ESO, I H 1 X . MAT . MAT2 . MATX . MX AX . 

1  NE  « NEAX . NEREM , NESC , NESX . NPREM , NR , NR 1  , NRAN , NR AX . NRREM , NRSO , 

2  NRSX,NRWL,NSCT,NSEG.NSPL,NTHX.NZ.NZ1 .  NZAX.NZBT. NZREM, NZSO 

3  NZSX.NZTP. PART, RAN, RHO, RH02. RSO, S. SPEC, SPEC2.SPH.SUZ.UC.UT 

4  UX , UXREM , UY , UYREM , UZ , UZREM , W , WMAX , WREM , WT , WTREM , X , XN , XREM . 

5  XSEC , XSECB , Y . YN , YREM , Z , ZH , ZN , ZREM , ZSO 


MATOR 


COMMON  CTHO, I WCT.NASX.SOANG 
IF (UZ ) 1 10,70* 10 
10  S8= (ZH(NZTP)— Z) /UZ 
ZNdZH (NZTP) 

XN«UX*SB+X 

YNfcUY*SB+Y 

R2=XN**2+YN**2 


NZSNZTP 

I F ( R2-RH02 ( NRWL ) ) 30 , 20 , 70 
20  NRcNRWL 
GO  TO  120 

30  DO  40  NRC=1.NRWL 
NR « NR C 


I F ( R  2-RH02 ( NRC ) >60,50,40 
40  CONTINUE 

50  IF (UX*XN+UY*YN )60 , 120  -  120 
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60  NR=NR-1 
GO  TO  120 
70  NR=NRWL 

DRADnUX*X+UY*r' 

SUZ2»UX**2+UY**2 

S8= ( SORT ( DRAD**2-SU22* ( X**2+Y**2-RH02 ( NRWL > ) } -DRAD ) /SUZ2 

XN=UX*SB+X 

YN=UY*S8+Y 

ZNaUZ*SB+Z 

DO  80  N2C  =  N2BT  .NZTp 

NZ“NZC 

IF  (2N-2H(N2C ) ) 100 » 90 » 80 
80  CONTINUE 
90  IF(U2) lOOi 120, 120 
100  N2=N2— 1 

GO  TO  120 

110  S8=(2H(N2BT)-2)/U2 
2N=2H(N2BT) 

XN=SB*UX+X 
YN*iSB*UY+Y 
R2=XN**2+YN**2 
N2=N2BT - 1 

I F ( R2-RH02 ( NRWL ) ) 30 . 20 . 70 
120  X=XN 
Y«YN 
2aZN 
NSGT  «  0 
RETURN 
END 


$  FASTRAN 

C  SUBROUTINE  C0MPT2  3-2-67 

C  SAMPLE  NEW  DIRECTION  AND  ENERGY  FROM  COMPTON  DISTRIBUTION 

SUBROUTINE  COMPT 2 

DIMENSION  CPU (360) .CTh(20>  4  DEN ( 4 )  ,EB(25> ,ECLS(30>  4  MAT ( 10.10)4 

1  MAT2( 10. 10) 4  NEREM ( 5o )  4  NPREM ( 50 ) 4  NRREM { 50 ) 4NSPL( IO4 10) 4 

2  N2REM ( 50 ) 4  PART (  10 4 10)  4 RHO ( 10) 4RH02( 10) .RSO(  10)  4  SPEC (20 4 30)  4 

3  SPE  C2  (  20  4  30  )  .SPH(360)  4UC000.4)  ,UT(300  44)  4UXREM(50)  4 

4  UYREM (50)4  U2REM ( 50  >  4  WREM (50)4  WTREM (50). XREM ( 50 ) 4  XSEC ( 2  4  4 ) 4 

5  XSECB ( 25  4  2  4  4 ) . YREM( 50  >  4 ZH ( 1 0 ) 4 ZREM<50) .ZSO(  10 ) 

DIMENSION  ESO<2) .CTh0( 10) 4S0ANG( 10)4 IWCT(4) .WMAX(4) 

COMMON  CPH  4  C.Th  4  DEN  4  EB  4  ECLS  4  ESO  4  I H I  X . MAT ,MAT2 4 MATX 4 MXAX 4 

1  NE  4  NEAX  4  NEREM, NESO , NES  X . NPREM , NR , NR  1  , NRAN , NR AX  4  NRReM , NRSO  4 

2  NRSX  4  NRWL  «  NSCT  4 NSEG4  NSPL . NTHX  4  NZ  4  NZ 1  4  NZAX 4 NZBT 4 NZREM , NZSO 4 

3  NZSX4N7TP.PART4RAN.RHO ,RH024RS0,S4SPEC4SPEC2.£PH,SUZ4UC4UT . 

4  UX  4  UXREM  4  UY . UYREM . UZ . UZREM , W  4  WMAX  4  WREM  4  WT  4  WTREM  4  X  4  XN . XREM  4 

5  XSEC  4  XSECB . Y  4  YN . YREM 4  Z  4  ZH . ZN  4  ZREM 4  ZSO 
COMMON  ctho. 1WCT.NASX,S0ANG 

10  RAN=RAM2B ( 0 ) 

T--2«  0/W 

I F ( RAN— ( l *  0  +  T ) / ( 9  4  0  +T )  ) 20 4 20 4 30 
20  RAN=R AM2B ( 0 ) 

R  =  1 «0  +  RAN*T 
RAN=RAM2B(0 ) 

IF (RAN-4 ,0*(R~1 .0)/(R*F2 ) >40.40. 10 


O  {3 
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30  RAN=RAM2B(0) 

Rs{ 1  * 0+T ) / ( 1  « 0+RAN*T > 

RAN=R AM28  <  0 ) 

IF  (  RAN— 0  »  5*  <  (  W-R*W+1  .  0 ) **2+1  . 0/R)  >40  <  '(0i  10 
40  WNs W*R 

C0M  = 1  .0+W-WN 
W=  WN 

NMAT=MAT2(NR.NZ) 

IF <W-WMAX(NMAT > >45.80,30 
45  SOM=SQRT ( 1 .0— C0M**2 ) 

NE= 1 0  2 • 2/W+O  «  5 
RAN=RAM2B(0) 

IPH=360.0*RAN 

UZN=UZ*COM+SUZ*SOM*CPH(  I  PH+1  ) 

SUZN=  SORT  (  1  »0-UZN**2) 

A=:SUZ*5UZN 

IFIA-0.000001 ) 50 , 50, 60 
50  UXN=- 3PH ( I PH+ 1 ) *SUZN 
UYN=SPH( IPH+1 )*SUZN 
GO  70 

60  CDPH= <COM-UZ*UZN)/A 

SDPHcSOM*SPH( IPH+1 )/SUZN 
UXNa( ( UX*CDPH— UY*SDPH ) *SUZN ) /SUZ 
UYN= (  <UY*CDPH+UX*oOPH) *SUZN) /SUZ 
70  UX=UXN 
UY»UYN 
UZ  n  UZN 
3UZ=SUZN 
30  RETURN 
END 


FASTRAN 

SUBROUTINE  GRADE  3-2-67 

C  RECORDS  SCORE  WHEN  PHOTON  ENTERS  DETECTOR 

SUBROUTINE  GRADE 

DIMENSION  CPHI360)  . CTh(20)  .  DEN ( 4 )  , ESC  25)  .ECLS(30>  « MAT ( 10.10). 

1  MAT2 (10.10) . NEREM ( 50 >  .NPREMC  50 )  . NRREM ( 50 ) »NSPL (10.10). 

2  NZREM ( 50 )  «  PART (  10.10)  . RHO ( 10) .RH02( 10) .RSO(  10) .SPEC (20. 30) . 

3  SPEC2<  20.30)  .SPH( 360 )  .UC ( 300.4  )  .UT (300.4 )  ,UXR£M( 50)  . 

4  UYREM(SO)  .UZREMC50)  . WREM ( 50 )  . WTREM(50)  .XREM(50) .XSEC(2.4  )  , 

5  XSEC8( 25.2.4) . YREM ( 50 ) . ZH ( 1 0 ) .ZREM(50) .ZS0( 10) 

DIMENSION  ESO ( 2 ) .CThOC 10) .SOANGl 10). IWCT(4 ) , WMAX( 4 ) 

COMMON  CPH . CTh . DEN . EB . ECLS . ESO , I H I  X . MAT , MAT2 . M ATX ♦ MXAX . 

1  NE . NEAX . NEREM . NESO ♦ NESX . NPREM . NR , NR  1  , NRAN . NR AX . NRREM , NRSO . 

2  NRSX.MRWL.NSCT.NSEG.NSPL.NTHX.NZ.NZl ,  NZAX.NZBT . NZREM, NZSO. 

3  NZSX , NZTP , PART , RAN , RHO , RH02 , RSO , S . SPEC . 3PEC2 . SPH . SUZ . UC ♦ UT . 

4  UX , UXREM . UY ♦ UYREM , UZ . UZREM , W . WMAX . WREM . WT . WTREM . X . XN ♦ XREM . 

5  XSEC . XSECB . Y , YN . YREM . Z . ZH . ZN , ZREM . ZSO 
COMMON  CTHO, IWCT.NASX.SOANG 

DO  10  NTH= 1 «NTHX 
JTHaNTH 

IF (UZ-CTH(NTH> ) 10.20.20 
10  CONTINUE 
20  E  =  0  «  5 1  1 XW 

DO  30  NEC=1.NEAX 


o  o 
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NECL=NEC 

IF  <  E-ECLS ( NEC )  >30*40,40 
30  CONTINUE 

40  SPEC  (  JTH , NECL ) =  SPEC ( JTH  *  NECL ) +WT 

SPEC 2 ( UTH  »  NECL )=SPEC2(UTH,  NECL ) +WT**2 

RETURN 

END 


$  FASTRAN 

SUBROUTINE  TEACH  4-20-67 

NORMALIZES  AND  PRINTS  RESULTS 
SUBROUTINE  TEACH 

D I  MENS  I  ON  CPH ( 360 )  * CTh ( 20 )  * DEN ( 4 )  , EB  < 25 ) , ECLS  (  30 )  * MAT (10*10)* 

1  MAT2< lOi 10) * NEREM ( 50 ) ,NPREM(50) ,NRREM(50) *NSPL( 10* 10) * 

2  NZREM ( 50 )  * PART  <  1 0  , 1 0 )  * RHO (10)*  RH02 (10)* RSO ( 1 0 ) «  SPEC (20*30)* 

3  SPEC2 ( 20  *  30 )  *SPH(360)  «UC(300*4)  *UT( 300*4) *UXREM(50) * 

4  UYREM ( 50 )  *  UZREM ( 50 )  *WREM(50)  ,WTREM(50)  * XREM ( 50  )  *XSEC(2*4 ) * 

5  XSECB ( 25  *2*4 ) «  YREM ( 50 )  * ZH ( 1 0 )  *ZREM(50)  *ZSO(  10) 

DIMENSION  ES0(2) *CTh0< 10) *SOANG( 10) . IWCT(4) »WMAX(4) 

COMMON  CPH * CTH* DEN, EB* ECLS, ESO, I H I  X , MAT , MAT2 * MATX » MXAX , 

1  NE . NEAX  *  NEREM , NESO , NESX  *  NPREM , NR , NR  1  , NRAN , NR AX  *  NRREM , NRSO , 

2  NRSX  *  NRWL  ,  NSCT  «  NSEG*  NS  PL  *IJTHX,NZ,NZ1  *  NZAX  *  NZBT  ,  NZREM,  NZSO  , 

3  NZSX  ,  NZTP  ,  PART  ,  RAN  ,  RHO  ,  RH02  ,  RSO  ,  S  ,  SPEC  ,  SPEC-2  ,  SPH  ,  SUZ  ,  UC  ,  UT  . 

4  UX , UXREM , UY , UYREM , UZ  *  UZREM , W , WMAX , WREM , WT , WTREM , X , XN , XREM , 

5  XSE  C , XSECB  *  Y , YN , YREM , Z . ZH  *  ZN  * ZREM , ZSO 
COMMON  CThO, I  WCT « NASX , SOANG , LP 

DIMENSION  ASPEC ( 20 ) ,ASPEC2(20) ,ESPEC(30) *ESpEC2(30) 

W0T6, 10 

10  FORMAT (  1  HO ) 

FLP=LP 
HI  X= I HIX 
M  I  X«=H  I  X*FLP 
DO  11  NEC=1,NEAX 
ESPEC (NEC) =0.0 
ESPEC2 ( NEC ) =0.0 
DO  11  NTH=1*NTHX 

ESPEC ( NEC ) cESPEC ( NEC ) +SPEC ( NTH  *NEC ) 

1 1  ESPEC2(NEC)=ESPEC2(NEC)+SPEC2(NTH*NEC) 

DO  1 2  NTH=  1  ,NTHX 
ASPEC (NTH) =0.0 
ASPEC2(NTH)=0.0 
DO  12  NEC=1,NEAX 

ASPEC ( NTH ) = ASPEC ( NTH ) +SPEC ( NTH  *  NEC ) 

12  ASPEC2(NTH) =ASPEC2(NTH)+SPEC2(NTH,NEC) 

DO  15  NTH=l*NTHX 
ASPEC (NTH) = ASPEC ( NTH) /HIX 
A SPEC 2 (NTH) =ASPEC2(NTH)/HIX 

ASPEC  2 ( NTH )  =  <  A SPEC 2 ( NTH ) -ASPEC ( NTH ) **2 ) / ( H IX- 1 .0) 

IF (ASPEC (NTH) )  13*  14 , 13 

13  ASPEC2 (NTH) =SQRT( ASPEC2(NTH) ) /ASPEC (NTH) 

GO  TO  15 

14  ASPEC2(NTH) =-10.0 

15  CONTINUE 
DO  18  NEC= 1 , NEAX 
ESPEC ( NEC ) =ESP£C ( NEC ) /H  I  X 
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ESPEC2<N£C)cESPEC2(NEC)/HIX 

E SPEC 2 ( NEC )  =  ( E SPEC 2 ( NEC ) -ESPEC  <  NEC )#*2)/(HIX-l  .0) 

I F ( ESPEC ( NEC ) ) 16, 17,16 

16  ESPEC2(NEC)=SQRT<ESPEC2(NEC) ) /ESPEC (NEC) 

GO  TO  18 

17  ESPEC 2 ( NEC ) =- 1 0 • 0 

18  CONTINUE 

DO  40  NTH= 1  *  NTHX 
DO  40  NEC  = 1  .NEAX 

SPEC  (  NT H  4 NEC  )  t=SPEC  ( NTH ,  NEC  )/H  1  X 
SPEC 2 (NTH*  NEC ) =SPEC2 (NTHi NEC ) /H I  X 

SPEC 2 ( NTH i NEC) =(SPEC2(NTH.NEC)-SPEC(NTH»NEC)**2)/<HIX-1 .0) 

I F  <  SPEC ( NTH , NEC )  ) 20 , 30  ,  20 

20  SPEC2 ( NTH , NEC ) =SQRT ( SPEC2 ( NTH  ,  NEC )  ) /SPEC ( NTH  *  NEC ) 

GO  TO  40 

30  SPEC2 ( NTH  *  NEC ) =- 1 0 • 0 
40  CONTINUE 
W0T6  4  50 

50  FORMAT (46H  ENERGY-ANGLE  SPECTRA  EMERGING  FROM  COLLIMATOR) 

W0T6  4  60 

60  FORMAT <65h  CTh  IS  THE  COSINE  OF  THE  LARGEST  ANGLE  IN  GIVEN  ANGULAR 
1  INTERVAL) 

W0T6  4  70 

70  FORMA"  :*»9H  ECLS  IS  SMALLEST  ENERGY  IN  GIVEN  ENERGY  INTERVAL) 

WOT61 80 

80  FORMAT  (46H  DATA  ARE  NORMALIZED  TO  ONE  PHOTON  ~RU,-,  SOURCE) 
a  0  T  6  4  1  0 

I F ( NEAX- 10)84,84482 
82  NBOT  = 1 
NT  0P=  1  0 
GO  TO  86 
84  NBO'f  =  1 

NT  0P=  NEAX 

86  WOT649O4 (ECLS(NEC) , NEC=NBCT ,NTOP ) 

90  FORMAT ( 1 2H  CTH,  ECLS=F7 . 3 , 9F 1 0 . 3 ) 

W0T6  4  1 0 

DO  110  NTH= I , NTHX 

W0T6,  100,CTH(NTH)  ,  ( SPEC ( NTH , NEC ) 4  NEC  =NBOT , NT OP ) 

100  FORMAT (F7. 3, 1  PE  1 4 . 2 , 1P9E10.2) 

1  1 0  CONT I NUE 
W0T6  4  1  0 

W0T6, 1 1 1 , ( ESPEC ( NEC ) , NEC=NQ0T , NT OP ) 

111  FORMAT ( 7H  SUM  1  PE  1 4 . 2 4  1 P9E 1 0 .2 ) 

W0TS4 1 0 

IF(NEAX-NTOP) 1164116,112 
1  12  NBOT  =  NTOP+ 1 
NT  OP  =  NTOP+ 1 0 
I F ( NEAX-NTOP  >114  4  86486 
114  NT  0P=  NEAX 
GO  TO  86 
1 16  WQT64 1 20 

120  FORMAT ( 59H  FRACTIONAL  STATISTICAL  DEVIATION  OF  SPECTRA  IN  ABOVE  TA 
1  BLE) 

WOT64 10 

IF ( NEAX- 1 0 )  124,124,122 
122  NBOT  =  1 
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NT  0P= 1 0 
GO  TO  126 
124  N30T  =  1 

NT  OP=  NEAX 

126  W0T6  *  90 *  ( ECLS ( NEC ) *  NEC  =NBOT *NTOP ) 

WO  T6  4  10 

DO  130  NTH=1*NTHX 

WO  T6t  140*CTH(NTH)  ,  < SPEC2 < NTH ♦ NEC )  , NEC= NBOT , NTOP > 
130  CONTINUE 

140  FORMAT (F7.3,Fl4.4*9F10.4  ) 

WOT 6  *  1 0 

W0T6  *  141  *  ( ESPEC2 ( NEC) , NEC  =  NBOT , NTOP ) 

141  FORMATtlOH  SUM  DEV.  FI  1 .4.9F10.4 > 

W  0  T  6  *  10 

IF < NEAX-NTOP) 146*146*142 

142  N80T=NT0P+1 
NT  0P=  NT  0P+ 1 0 

IF ( NEAX-NTOP ) 144* 126*126 
144  NT  0P  =  NEAX 
GO  TO  126 
146  W0T6*150 

150  FORMAT ( 30H  ANGULAR  FRACTIONAL) 

W0T6  *  1 60 

160  FORMA  T ( 29H  CTH  DISTRIBUTION  DEVIATION) 

DO  170  NTH= 1  *  NTHX 

WOT 6  *  1 80  *  CT H ( NTH )  * ASPEC(NTH)  *ASPEC2<NTh) 

170  CONTINUE 

180  FORMAT (F7. 3, 1PE12.2,0PF1 1.4) 

WOT 6 «  10 
RE  TURN 
END 


S  FASTRAN 

C  SUBROUT  INI'.  TAB  IN  28-3-64 

SU8R0UT INE  T AB I N ( NT AB I N * FB * XB * NMAX * MMAX *  X  * FX ) 

DIMENSION  FB(25*8) *X8(25) *FX(Q) *XBAV(25) *D1 <25) ,D2(25) 
GO  10  (10*30) *NTABIN 
10  NMAX 1 =NMAX— 1 

DO  20  N=2  *  NMAX 1 

X8A V ( N ) = ( XB ( N— 1 )+XS(N) )/2*0 

D 1 ( N ) = ( XB ( N- 1 ) -XB ( N ) ) * ( XB ( N- 1 ) -XB ( N+ ) > ) 

20  D2(N>  =<XB(N)-XB(N— 1  )  ) * ( XB ( N ) -XB < N+ 1  )  ) 

NTAB I N=2 
NEXSsNMAXl-2 
30  IF (X-X8( 1 ) >60,50*40 
40  NX  =2 

GO  TO  200 
50  NX  = 1 

GO  TO  220 

60  IF ( X-XB( 2 ) ) 80,70* 40 
70  NX  =2 

GO  TO  220 

80  IF(NEXS)90* 110*140 
90  W0T6* 1 00 

100  FORMAT (33HONOT  ENOUGH  BASE  POINTS  FOR  T AB I N ) 
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CALL  SYSTEM 

110  IF(X-XB(NMAX) > 130. 120, 130 
120  MX  =  NM A X 

GO  TO  220 
130  NX  =  NM  AX  1 
GO  TO  200 

1  A 0  DO  170  N  =  3 ♦ NMAX 1 

I F ( X^XB ( N) ) 170*150.160 
150  NX=N 

GO  TO  220 
160  NX=N 

GO  TO  180 
170  CONTINUE 
GO  TO  110 

180  IF ( X-XBAVC NX )  >200,200, 190 
1O0  NX=NX-1 

200  'tl  T  1  =  (  X-XB(NX)  )  *  (  X— XQ  (  NX+  1  )  ) /D  1  (NX) 

W T 2= ( X-X3 ( NX-1 ) ) * ( X-XB ( NX+ 1 ) ) /D2 ( NX ) 

WT  3= 1  .0— WT1-WT2 
DO  210  M= 1 , MMAX 

210  F X ( M ) =  WT 1 *FB ( NX— 1  , M ) +WT2*F8 < NX , M ) +  WT 3*FB ( NX+ 1 ,M) 
RETURN 

220  DO  230  M  = 1  , MMAX 
230  F X ( M ) =FB  (  NX  «  M ) 

RETURN 

END 
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